!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module fv_restart_mod use fms_mod, only: file_exist, error_mesg, lowercase, FATAL, & NOTE, read_data, uppercase, stdout, mpp_pe, mpp_root_pe, mpp_chksum use fms_io_mod, only: open_restart_file, close_file, write_data, get_instance_filename use mpp_io_mod, only: mpp_open, mpp_close, axistype, fieldtype, & mpp_read_meta, mpp_get_info, mpp_get_fields, mpp_read, & MPP_NETCDF, MPP_SINGLE, MPP_MULTI, MPP_SEQUENTIAL, MPP_RDONLY, & MPP_WRONLY, MPP_NATIVE, MPP_OVERWR, mpp_write_meta, mpp_write, & atttype, mpp_get_atts, mpp_get_att_name, mpp_get_att_real_scalar, & mpp_get_field_name use fv_pack, only: ak, bk, ks, ptop, beglat, endlat, & master, nlon, mlat, nt_phys, nqrst, ng_d, ng_s, & ighost, mountain, coslon, sinlon, lon, lat, lonb, latb, & adjust_dry_mass, use_set_eta, p_var, d2a3d, & age_tracer, age_time, get_eta_level, & do_fms_tracer_manager, fms_tracers_file, fv_domain, & restart_format, use_tendency, cold_start, gid #ifdef MARS_GCM use fv_pack, only: p_ref #endif MARS_GCM use tracer_manager_mod, only : tr_get_tracer_names=>get_tracer_names, & get_tracer_names, get_number_tracers, & set_tracer_profile, & get_tracer_index, NO_TRACER use field_manager_mod, only : MODEL_ATMOS use platform_mod, only: I8_KIND ! use shr_kind_mod, only : r8 => shr_kind_r8 #ifdef SPMD use mod_comm, only: mp_bcst_int, mp_bcst_r2d, mp_scatter3d, mp_bcst_real, & mp_send4d_ns, mp_recv4d_ns, mp_gather3d #endif use fv_arrays_mod, only: fv_array_check, fv_stack_push, fv_array_sync, & fv_print_chksums, ksp, kep use mpp_mod, only: mpp_sync, mpp_root_pe, mpp_broadcast use mpp_domains_mod, only: mpp_get_data_domain, mpp_get_compute_domain use constants_mod, only: grav, kappa, rdgas use pmaxmin_mod, only: pmaxming, pmaxmin implicit none private public :: fv_restart, read_fv_rst, write_fv_rst, add_tracers contains subroutine read_fv_rst(im, jm, km, nq, dflnm, days, seconds, oform) #include "fv_arrays.h" character*(*), intent(in) :: dflnm ! dynamics filename integer, intent(in) :: nq integer, intent(out):: im, jm, km integer, intent(out):: days integer, intent(out):: seconds character(len=*), intent(in) :: oform ! Local variables integer iuic ! Unit number integer i, j, k, n, jn, it, ntracers logical is_tracer_in_restart_file character(len=128) :: tname real zsurf(nlon, beglat: endlat) real ginv, qmax, qmin real, allocatable :: w2d(:,:), w3d(:,:,:) type(fieldtype), allocatable :: tracer_fields(:), fields(:) type(atttype), allocatable :: attributes(:) integer :: ndim,natt,nvar,ntime real p0 character(len=*), parameter :: fms_tracers_restart_file = & 'INPUT/'//fms_tracers_file character(len=128) :: directory, restart_name, tr_name, dflnm_mod logical :: use_native_format integer :: ids, ide, jds, jde, ics,ice,jcs,jce #ifdef MARS_GCM ! ------Dont require specific water vapor fields in netcdf version ! of tracer restart file: Hence this code is not necessary #else integer :: nsphum, nliq_wat, nice_wat, ncld_amt character(len=20) :: fv_tr_name(4) data fv_tr_name / "sphum","liq_wat","ice_wat","cld_amt" / #endif #ifdef use_shared_pointers pointer( p_zsurf, zsurf ) #include "fv_point.inc" call fv_stack_push( p_zsurf, nlon*(endlat-beglat+1) ) #endif ginv = 1./ GRAV !Modify the filename if it has an ensemble instance appendix !The following call does not modify the filename unless there !is an ensemble instance appendix and fms_io knows about it. call get_instance_filename(dflnm, dflnm_mod) ! output format if (trim(uppercase(oform)) == 'NATIVE') then use_native_format = .true. restart_name = trim(dflnm_mod) else if (trim(uppercase(oform)) == 'NETCDF') then use_native_format = .false. if (trim(dflnm_mod(LEN_TRIM(dflnm_mod)-2:LEN_TRIM(dflnm_mod))) /= ".nc" ) & restart_name = trim(dflnm_mod)//".nc" if (file_exist(dflnm_mod) .and. .not. file_exist(trim(dflnm_mod)//".nc")) then call error_mesg( 'read_fv_rst', & 'NetCDF file does not exist. Switching to Native format.', & NOTE) restart_name = trim(dflnm_mod) use_native_format = .true. endif else call error_mesg ('read_fv_rst', & 'invalid value for argument oform='//trim(oform), FATAL) endif if (.not. file_exist(dflnm_mod) .and. .not. file_exist(trim(dflnm_mod)//".nc")) & call error_mesg ('read_fv_rst', & 'Neither a Native or NetCDF restart file is present.', FATAL) if( master ) then if (use_native_format) then if (file_exist(dflnm_mod)) then !!$ call mpp_open (iuic, dflnm_mod, action=MPP_RDONLY, form=MPP_NATIVE, access=MPP_SEQUENTIAL) iuic = open_restart_file(dflnm_mod, 'read') read(iuic) im, jm, km, nqrst endif else if (file_exist(trim(dflnm_mod)//'.nc')) then call mpp_open(iuic, dflnm_mod, & action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_SINGLE ) call mpp_read_meta(iuic) call mpp_get_info(iuic,ndim,nvar,natt,ntime) allocate(attributes(natt)) call mpp_get_atts(iuic,attributes) do n= 1,natt ! select case (attributes(n)%name) select case ( mpp_get_att_name(attributes(n))) case ('nlon') im= int(mpp_get_att_real_scalar(attributes(n))) case ('mlat') jm= int(mpp_get_att_real_scalar(attributes(n))) case ('nlev') km= int(mpp_get_att_real_scalar(attributes(n))) case ('nt_phys') nqrst= int(mpp_get_att_real_scalar(attributes(n))) case ('days') days= int(mpp_get_att_real_scalar(attributes(n))) case ('seconds') seconds= int(mpp_get_att_real_scalar(attributes(n))) case ('ks') ks = int(mpp_get_att_real_scalar(attributes(n))) end select enddo deallocate(attributes) allocate(fields(nvar)) call mpp_get_fields(iuic, fields) do n=1,nvar ! select case (fields(n)%name) select case ( mpp_get_field_name(fields(n))) case ('ak') call mpp_read(iuic,fields(n),ak) case ('bk') call mpp_read(iuic,fields(n),bk) end select enddo call mpp_close(iuic) deallocate(fields) endif endif if ( im /= nlon .or. jm /= mlat .or. km /= nlev ) then write(6,*) 'Dataset resolution:',im, jm, km, nqrst call error_mesg('read_fv_rst:','resolution inconsistent', FATAL) endif if ( nqrst > nt_phys .or. nqrst > nq) then write(6,*) ' ' write(6,*) '---- Warning -------' write(6,*) 'Total number of tracers in restart file is greater than either' write(6,*) 'NT_PHYS or NQ. The number NT_PHYS from the namelist is used.' write(6,*) 'NQRST=', nqrst, 'NT_PHYS=', nt_phys, ' NQ=', nq nqrst = nt_phys endif if (use_native_format) then read(iuic) days, seconds read(iuic) ak, bk, ks endif write(6,*) 'Restart days=', days,' sec=',seconds write(6,*) '----------------------------------' write(6,*) 'Model Top Pressure (pa)= ', ak(1) write(6,*) '----------------------------------' write(6,*) 'Checking Initial condition ...' endif #ifdef SPMD ! call mp_bcst_int(days) ! call mp_bcst_int(seconds) ! call mp_bcst_int(im) ! call mp_bcst_int(jm) ! call mp_bcst_int(km) ! call mp_bcst_int(nqrst) ! call mp_bcst_int(ks) ! ! call mp_bcst_r2d(km+1, 1, 1, 1, ak, 0) ! call mp_bcst_r2d(km+1, 1, 1, 1, bk, 0) call mpp_broadcast( days, mpp_root_pe() ) call mpp_broadcast( seconds, mpp_root_pe() ) call mpp_broadcast( im, mpp_root_pe() ) call mpp_broadcast( jm, mpp_root_pe() ) call mpp_broadcast( km, mpp_root_pe() ) call mpp_broadcast( nqrst, mpp_root_pe() ) call mpp_broadcast( ks, mpp_root_pe() ) call mpp_broadcast( ak, km+1, mpp_root_pe() ) call mpp_broadcast( bk, km+1, mpp_root_pe() ) #endif ptop = ak(1) if (use_native_format) then !--------------------------- ! read surface geopotential !--------------------------- allocate ( w2d(im,jm) ) #ifdef SPMD if ( master ) read(iuic) w2d call mp_scatter3d(w2d, phis, im, jm, 1, beglat, endlat, & 1, 1, ighost, ighost, 0) call fv_array_sync #else read(iuic) w2d do j=beglat,endlat do i=1,im phis(i,j) = w2d(i,j) enddo enddo #endif ! call pmaxming('PHIS', phis, im, jm, 1, beglat, endlat, & ! ighost, ighost, ginv) allocate ( w3d(im, jm, km) ) !--------------------------- ! Read u-wind !--------------------------- #ifdef SPMD if ( master ) read(iuic) w3d call mp_scatter3d(w3d, u, im, jm, km, beglat, endlat, & 1, km, ng_d, ng_s, 0) call fv_array_sync #else read(iuic) u #endif !--------------------------- ! Read v-wind !--------------------------- #ifdef SPMD if ( master ) read(iuic) w3d call mp_scatter3d(w3d, v, im, jm, km, beglat, endlat, & 1, km, ng_d, ng_d, 0) call fv_array_sync #else read(iuic) v #endif !----------------------------------- ! Read virtual potential temperature !----------------------------------- #ifdef SPMD if ( master ) read(iuic) w3d call mp_scatter3d(w3d, pt, im, jm, km, beglat, endlat, & 1, km, ng_d, ng_d, 0) call fv_array_sync #else read(iuic) pt #endif !----------------------------------- ! Read pressure thickness (pascal) !----------------------------------- #ifdef SPMD if ( master ) read(iuic) w3d call mp_scatter3d(w3d, delp, im, jm, km, beglat, endlat, & 1, km, 0, 0, 0) call fv_array_sync #else read(iuic) delp #endif else !this code may not work? all these fields do not have the same data domain call mpp_get_data_domain( fv_domain, ids, ide, jds, jde ) call mpp_get_compute_domain( fv_domain, ics, ice, jcs, jce ) call read_data ( restart_name, 'Surface_geopotential', phis(ids:ide,jcs:jce), domain = fv_domain) call read_data ( restart_name, 'U', u(ids:ide,jds:jde,:) , domain = fv_domain) call read_data ( restart_name, 'V', v(ids:ide,jds:jde,:), domain = fv_domain) call read_data ( restart_name, 'T', pt(ids:ide,jds:jde,:), domain = fv_domain) call read_data ( restart_name, 'DELP', delp(ids:ide,jds:jde,:), domain = fv_domain) endif ! not use_native_format #ifdef SHAVE_P call set_eta(km, ks, ptop, ak, bk) do j=beglat,endlat do i=1,im delp(i,j,1) = ak(2) - ak(1) enddo enddo #endif !------------------------------------ ! Read tracers one at a time !------------------------------------ if (use_native_format) then do n=1, nqrst #ifdef SPMD if ( master ) read(iuic) w3d call mp_scatter3d(w3d, q(1,beglat-ng_d,1,n), im, jm, km, beglat, endlat, & 1, km, ng_d, ng_d , 0) call fv_array_sync #else read(iuic) w3d do k=1,nlev do j=beglat,endlat do i=1,im q(i,j,k,n) = w3d(i,j,k) enddo enddo enddo #endif enddo else ! ----------read "dynamics" tracers from tracer restart file #ifdef MARS_GCM !rjw ------No longer require specific water vapor fields in netcdf version of tracer restart file !rjw It now makes more sense to read all the tracers from the tracer file !rjw in a single loop, that appears below in the processing of field manager tracers #else do n =1,nqrst call read_data(fms_tracers_restart_file,fv_tr_name(n),q(ids:ide,jds:jde,:,n) , domain = fv_domain) enddo #endif MARS_GCM endif ! --------------- Finished reading dynamcics tracers ------------ if( age_tracer ) call error_mesg('read_fv_rst', 'age_tracer == .T. not yet working', FATAL) if (use_native_format) then if ( age_tracer .and. nqrst == nq ) then ! this check will not work properly (pletzer) ! q = 0. ; age_time = 0. if(master) read(iuic) age_time #ifdef SPMD ! call mp_bcst_real(age_time) call mpp_broadcast( age_time, mpp_root_pe() ) #endif endif if ( master ) call close_file(iuic) endif !----------------------------------------- ! Read additional tracers from netCDF file !----------------------------------------- #ifdef MARS_GCM if (do_fms_tracer_manager ) then !{ #else if (do_fms_tracer_manager .and. nq > nt_phys) then !{ #endif MARS_GCM !Modify the filename if it has an ensemble instance appendix call get_instance_filename(fms_tracers_restart_file,restart_name ) restart_name=trim(restart_name) if( .not. file_exist(restart_name//'.nc') ) then !{ if(master) then write(*,*)'---- Warning -------' write(*,*)'do_fms_tracer_manager=',do_fms_tracer_manager, & ' no of tracers nq=', nq write(*,*)'but netcdf restart file ', restart_name//'.nc', ' does not exist' write(*,*)'--------------------' endif else !}{ if(master) then write(*,*)'Additional tracers present in field_table. ' write(*,*)'These will be read from netcdf file ',restart_name endif call mpp_open(iuic, restart_name, & action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_SINGLE ) call mpp_read_meta(iuic) call mpp_get_info(iuic,ndim,nvar,natt,ntime) if ( nt_phys+nvar /= nq ) then if (master) then write(6,*) ' ' write(6,*) '---- Warning -------' write(6,*) 'Total number of tracers nq from netcdf and binary restart files does not match' write(6,*) 'nt_phys=', nt_phys,' nvar=', nvar,' nq=', nq, ' nt_phys+nvar /= nq' write(*,*)'---------------------' endif endif allocate(tracer_fields(nvar)) !if(master) call mpp_get_fields(iuic, tracer_fields) ! read only the tracers that are registered in the tracer manager call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) n = 0 call mpp_get_data_domain( fv_domain, ids, ide, jds, jde ) do it = 1, ntracers !{ ! loop over tracers call get_tracer_names(MODEL_ATMOS, it, tname) #ifdef MARS_GCM #else if (use_native_format) then ! exclude tracers read from native format file if(tname == "sphum" ) cycle if(tname == "liq_wat") cycle if(tname == "ice_wat") cycle if(tname == "cld_amt") cycle endif #endif MARS_GCM ! this tracer will be added n = n + 1 is_tracer_in_restart_file = .FALSE. do jn = 1, nvar !{ ! loop over netcdf file tracers ! if( lowercase( trim(tracer_fields(jn)%name) ) == lowercase( trim(tname) ) ) then !{ if( lowercase( trim(mpp_get_field_name(tracer_fields(jn))) ) == lowercase( trim(tname) ) ) then !{ if(master) print *,'tracer ', trim(tname), ' found!' ! ok found the tracer in the restart file is_tracer_in_restart_file = .TRUE. ! read data and store in array w3d ! call mpp_read(iuic, tracer_fields(jn), w3d) call read_data( restart_name, tname, & q(ids:ide,jds:jde,:,it), domain = fv_domain) endif !} if(is_tracer_in_restart_file) exit enddo !} if(.NOT. is_tracer_in_restart_file) then !{ ! missing tracer from restart file if(master) then write(6,*) '---- Warning -------' write(6,*) 'Tracer ', tname, ' was defined in field table but could not be found in restart file.' write(*,*)'---------------------' endif ! set default profile call set_tracer_profile (MODEL_ATMOS, it, & q(ids:ide,jds:jde,:,it)) endif !} enddo !} call mpp_close(iuic) deallocate(tracer_fields) endif !} endif !} if (use_native_format) deallocate ( w3d ) !.............................................................................................. do j=beglat,endlat do i=1,im zsurf(i,j) = phis(i,j) enddo enddo call pmaxmin('ZS', zsurf, qmin, qmax, im*(endlat-beglat+1), 1, ginv) if ( qmax > 1. ) then mountain = .true. else mountain = .false. endif call pmaxming('U ', u, im, jm, km, beglat, endlat, ng_d, ng_s, 1.) call pmaxming('V ', v, im, jm, km, beglat, endlat, ng_d, ng_d, 1.) ! call pmaxming('DELP ', delp, im, jm, km, beglat, endlat, 0, 0, 1.) call pmaxming('T ', pt, im, jm, km, beglat, endlat, ng_d, ng_d, 1.) do n=1, nqrst call pmaxming('Q ', q(1,beglat-ng_d,1,n), im, jm, km, beglat, endlat, ng_d, ng_d, 1.) enddo !----------------------- ! Add additional tracers: !----------------------- #ifdef SET_Q_TOP ! Set water vapor at top layer to climatology value (for stratosphere) do j=beglat, endlat do i=1,im q(i,j,1,1) = 3.E-6 enddo enddo #endif if( use_set_eta ) call set_eta(km, ks, ptop, ak, bk) call p_var(nlon, mlat, nlev, beglat, endlat, ptop, delp, ps, & pe, peln, pk, pkz, kappa, q, ng_d, ncnst, adjust_dry_mass ) !--------------------------------------------------- ! Read surface pressure & winds for coupling purpose !--------------------------------------------------- !Modify the filename if it has an ensemble instance appendix call get_instance_filename('INPUT/fv_srf_wnd.res',restart_name ) restart_name=trim(restart_name) if (file_exist(restart_name) .or. file_exist(trim(restart_name)//'.nc')) then if (use_native_format) then if( master ) then !!$ call mpp_open (iuic, restart_name, action=MPP_RDONLY, & !!$ form=MPP_NATIVE, access=MPP_SEQUENTIAL) iuic = open_restart_file(restart_name, 'read') read(iuic) im, jm if ( im /= nlon .or. jm /= mlat ) then write(6,*) 'Dataset resolution:',im, jm call error_mesg('read_fv_rst:','resolution inconsistent in fv_srf_wnd', FATAL) endif endif ! Read ps_bp (surface pressure before physics update): if(master) read(iuic) w2d #ifdef SPMD call mp_scatter3d(w2d, ps_bp, im, jm, 1, beglat, endlat, & 1, 1, 0, 0, 0) call fv_array_sync #endif ! Read u-comp if(master) read(iuic) w2d #ifdef SPMD call mp_scatter3d(w2d, u_srf, im, jm, 1, beglat, endlat, & 1, 1, 0, 0, 0) call fv_array_sync #endif ! Read v-comp if(master) read(iuic) w2d #ifdef SPMD call mp_scatter3d(w2d, v_srf, im, jm, 1, beglat, endlat, & 1, 1, 0, 0, 0) call fv_array_sync #endif if(master) call close_file(iuic) else ! not native format call mpp_open(iuic, restart_name, & action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_SINGLE ) if (master ) then call mpp_read_meta(iuic) call mpp_get_info(iuic,ndim,nvar,natt,ntime) allocate(attributes(natt)) call mpp_get_atts(iuic,attributes) do n= 1,natt ! select case (attributes(n)%name) select case ( mpp_get_att_name(attributes(n))) case ('nlon') im= int(mpp_get_att_real_scalar(attributes(n))) case ('mlat') jm= int(mpp_get_att_real_scalar(attributes(n))) end select enddo deallocate(attributes) endif call mpp_close(iuic) if ( im /= nlon .or. jm /= mlat ) then write(6,*) 'Dataset resolution:',im, jm call error_mesg('read_fv_rst:','resolution inconsistent in fv_srf_wnd', FATAL) endif call read_data ( restart_name, 'ps_bp', ps_bp(ids:ide,jcs:jce), domain = fv_domain) call read_data ( restart_name, 'U_srf', u_srf(ids:ide,jcs:jce), domain = fv_domain) call read_data ( restart_name, 'V_srf', v_srf(ids:ide,jcs:jce), domain = fv_domain) endif else if(master) then write(*,*) ' ' write(*,*) '--------------------------------------------' write(*,*) 'File not found: ',restart_name write(*,*) 'Surface winds will be generated from D grid' write(*,*) '--------------------------------------------' endif do j=beglat,endlat do i=1,im ps_bp(i,j) = ps(i,j) enddo enddo call d2a3d(u(1,beglat-ng_d,km), v(1,beglat-ng_d,km), u_srf, v_srf, & nlon, mlat, 1, beglat, endlat, ng_d, ng_s, coslon, sinlon) endif call pmaxming('u_srf ', u_srf, im, jm, 1, beglat, endlat, 0, 0, 1.) call pmaxming('v_srf ', v_srf, im, jm, 1, beglat, endlat, 0, 0, 1.) if (use_native_format) deallocate ( w2d ) ! check sum input data call print_check_sum ('Check sums for FV input data:') if( master ) then call error_mesg ('read_fv_rst', 'FV netcdf restart file read.', NOTE) endif #ifdef PLET ! Do Pure Lagrangian Error Test if ( nq >= 7 ) then p0 = (1.E5) ** kappa !$omp parallel do private(i,j,k) do k=ksp,kep do j=beglat, endlat do i=1,im q(i,j,k,nq-2) = pt(i,j,k) q(i,j,k,nq-1) = pt(i,j,k)/pkz(i,j,k) * p0 q(i,j,k,nq) = pkz(i,j,k) * p0 enddo enddo enddo call fv_array_sync() endif #endif end subroutine read_fv_rst subroutine add_tracers( q, im, jm, km, nq, nqrst, beglat, endlat, ng, & age_tracer, age_time ) integer, intent(in):: im,jm,km,nq integer, intent(in):: nqrst integer, intent(in):: beglat, endlat integer, intent(in):: ng logical, intent(in):: age_tracer real, intent(inout):: age_time real, intent(inout):: q(im, beglat-ng:endlat+ng, km, nq) ! Local: integer i,j,k,iq call fv_array_check( LOC(q) ) do iq=nqrst+1, nq !$omp parallel do private(i,j,k) do k=ksp,kep do j=beglat, endlat do i=1,im q(i,j,k,iq) = 0. enddo enddo enddo enddo call fv_array_sync() if ( age_tracer ) age_time = 0. end subroutine add_tracers subroutine fv_restart(days, seconds) #include "fv_arrays.h" integer, intent(out):: days integer, intent(out):: seconds real, allocatable :: g_phis(:,:) real, allocatable :: qtmp(:,:,:) integer :: iuhs ! Unit to read in terrain data for cold start integer :: i,j,k integer :: im, jm, km character(len=128) :: restart_name #ifdef MARS_GCM character (len=128) :: filename, fieldname #endif MARS_GCM #include "fv_point.inc" call fv_print_chksums( 'Entering fv_restart' ) if (cold_start) then #ifdef SW_DYN if(nlev /= 1) call error_mesg('FV_init:','Set nlev=1 for Shallow Water Dynamics', FATAL) !----------------------- ! Shallow water dynamics !----------------------- ptop = 0. ak(1) = 0. ak(2) = 0. bk(1) = 0. bk(2) = 1. KS = 0 #ifndef USE_LIMA n_spong = 0 #endif call init_sw_ic(icase, nlon, mlat, nlev, ncnst, beglat, endlat, ng_d, ng_s, & phis, u, v, pt, delp, q, grav, radius, omega, f_d, ighost ) #else ! Start from a "dry & cold" isothermal atmosphere at rest #ifdef MARS_GCM filename= 'INPUT/surf.res.nc' fieldname= 'topo' if (file_exist( trim(filename) )) then mountain = .true. if(master) then allocate ( g_phis(nlon,mlat) ) call read_data( trim(filename), trim(fieldname), g_phis, no_domain=.true. ) endif call mp_scatter3d(g_phis, phis, nlon, mlat, 1, & beglat, endlat, 1, 1, ighost, ighost, 0) call fv_array_sync if(master) deallocate ( g_phis ) phis = phis * grav else mountain = .false. if(master) write(6,*) 'No surface dataset found; setting phis = 0' phis = 0. endif #else !Modify the filename if it has an ensemble instance appendix call get_instance_filename('INPUT/surf.res',restart_name ) restart_name=trim(restart_name) if (file_exist(restart_name)) then mountain = .true. if(master) call mpp_open(iuhs, restart_name, action=MPP_RDONLY, & form=MPP_NATIVE, access=MPP_SEQUENTIAL) #ifdef SPMD if(master) then allocate ( g_phis(nlon,mlat) ) read(iuhs) g_phis endif call mp_scatter3d(g_phis, phis, nlon, mlat, 1, & beglat, endlat, 1, 1, ighost, ighost, 0) call fv_array_sync if(master) deallocate ( g_phis ) #else read(iuhs) phis #endif if(master) close (iuhs) phis = phis * grav else mountain = .false. if(master) write(6,*) 'No surface dataset found; setting phis = 0' phis = 0. endif #endif MARS_GCM call set_eta(nlev, ks, ptop, ak, bk) call init_dry_atm(mountain, kappa, grav, rdgas) if ( age_tracer ) then age_time = 0. !$omp parallel do private (i, j, k) ! do k=1,nlev do k=ksp,kep do j=beglat, endlat do i=1,nlon q(i,j,k,ncnst) = 0. enddo enddo enddo call fv_array_sync() endif #endif else !--------------------------- ! Read from the restart file !--------------------------- restart_name ='INPUT/fv_rst.res' !The call get_instance_filename() is inside the read_fv_rst subroutine call read_fv_rst(im, jm, km, ncnst, restart_name, days, seconds, restart_format) ! call p_var(nlon, mlat, nlev, beglat, endlat, ptop, delp, ps, & ! pe, peln, pk, pkz, kappa, q, ng_d, ncnst, adjust_dry_mass ) endif #ifdef SPMD ! Ghost phis boundary N/S if ( ighost /= 0 ) then call mp_send4d_ns(nlon, mlat, 1, 1, beglat, endlat, & 1, 1, ighost, ighost, phis) call mp_recv4d_ns(nlon, mlat, 1, 1, beglat, endlat, & 1, 1, ighost, ighost, phis) call fv_array_sync() endif #endif !------------------------------------------------ ! Safety check of the (a,b) of the eta coordinate !------------------------------------------------ call check_eta(nlev, ks, ptop, ak, bk) ! Generate winds on the mass grid for the physics/coupler: if ( use_tendency ) then call d2a3d(u, v, ua, va, nlon, mlat, nlev, beglat, endlat, & ng_d, ng_s, coslon, sinlon) endif call fv_print_chksums( 'Exiting fv_restart' ) end subroutine fv_restart ! $Id: fv_restart.F90,v 17.0 2009/07/21 02:53:25 fms Exp $ subroutine write_fv_rst( dflnm, days, seconds, grav, oform) #include "fv_arrays.h" ! !INPUT PARAMETERS: integer, intent(in):: days integer, intent(in):: seconds real, intent(in):: grav character*(*), intent(in):: dflnm ! dynamics restart filename character(len=*), intent(in) :: oform ! local integer :: iuic ! Unit number real, allocatable :: w2d(:,:), w3d(:,:,:), latu(:), phis_alloc(:,:) real :: ginv, qmax, qmin real :: pfull(nlev), phalf(nlev+1) integer i, j, k, n character(len=128) :: tracer_name, tracer_longname, tracer_units type(axistype) :: axis_lon, axis_lat, axis_lev, axis_levh, & axis_lonv, axis_latu, axis_time type(fieldtype), allocatable :: tracer_fields(:) type(fieldtype) :: u_field, v_field, t_field, phis_field, delp_field type(fieldtype) :: ak_field, bk_field character(len=*), parameter :: fms_tracers_restart_file = & 'RESTART/'//fms_tracers_file character(len=128) :: directory, restart_name, tr_name logical :: use_native_format integer :: ids,ide, jds,jde, ics,ice,jcs,jce integer :: nsphum, nliq_wat, nice_wat, ncld_amt, ntracers real :: dummy(1) =1.0 #ifdef MARS_GCM !rjw fv_tr_name is not used in this subroutine #else character(len=20) :: fv_tr_name(4) data fv_tr_name / "sphum","liq_wat","ice_wat","cld_amt" / #endif MARS_GCM #include "fv_point.inc" ! output format if (trim(uppercase(oform)) == 'NATIVE') then use_native_format = .true. else if (trim(uppercase(oform)) == 'NETCDF') then use_native_format = .false. else call error_mesg ('write_fv_rst', & 'invalid value for argument oform='//trim(oform), FATAL) endif ! check sum output data call print_check_sum ('Check sums for FV output data:') if (use_native_format) then if ( master ) then !!$ if (file_exist(dflnm)) then !!$ call mpp_open ( iuic, dflnm, action=MPP_OVERWR, & !!$ form=MPP_NATIVE, access=MPP_SEQUENTIAL ) !!$ else !!$ call mpp_open ( iuic, dflnm, action=MPP_WRONLY, & !!$ form=MPP_NATIVE, access=MPP_SEQUENTIAL ) !!$ !!$ endif !Modify the filename if it has an ensemble instance appendix call get_instance_filename(dflnm,restart_name ) restart_name=trim(restart_name) iuic = open_restart_file(restart_name, 'write') write(iuic) nlon, mlat, nlev, nt_phys write(iuic) days, seconds write(iuic) ak, bk, ks write(6,*) ' ' write(6,*) 'FV restart file header written for days=', days, 'sec=',seconds endif call pmaxmin('PS', ps(1,beglat), qmin, qmax, nlon, endlat-beglat+1, 0.01) !--------------------------- ! Write surface geopotential !--------------------------- #ifdef SPMD allocate ( w2d(nlon,mlat) ) call mp_gather3d(phis, w2d, nlon, mlat, 1, beglat, endlat, 1, 1, & ighost, ighost, 0) if ( master ) then write(iuic) w2d ginv = 1./ grav call maxmin_global('Surface height', w2d, nlon, 1, mlat, ginv) endif #else write(iuic) phis #endif allocate ( w3d(nlon,mlat,nlev) ) !--------------------------- ! Write u-wind !--------------------------- #ifdef SPMD call mp_gather3d(u, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, ng_d, ng_s, 0) if ( master ) then write(iuic) w3d call maxmin_global('U', w3d, nlon, mlat, nlev, 1.) endif #else write(iuic) u #endif !--------------------------- ! Write v-wind !--------------------------- #ifdef SPMD call mp_gather3d(v, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, ng_d, ng_d, 0) if ( master ) then write(iuic) w3d call maxmin_global('V', w3d, nlon, mlat, nlev, 1.) endif #else write(iuic) v #endif !------------------------------------ ! Write virtual potential temperature !------------------------------------ #ifdef SPMD call mp_gather3d(pt, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, ng_d, ng_d, 0) if ( master ) then write(iuic) w3d call maxmin_global('T', w3d, nlon, mlat, nlev, 1.) endif #else write(iuic) pt #endif !------------------------------------ ! Write pressure thickness (pascal) !------------------------------------ #ifdef SPMD call mp_gather3d(delp, w3d, nlon, mlat, nlev, beglat, endlat, 1, nlev, 0, 0, 0) if ( master ) then write(iuic) w3d ! call maxmin_global('DELP', w3d, nlon, mlat, nlev, 1.) endif #else write(iuic) delp #endif if(ncnst /= 0) then !------------------------------------ ! Write tracers one at a time !------------------------------------ do n=1,nt_phys #ifdef SPMD call mp_gather3d(q(1,beglat-ng_d,1,n), w3d, nlon, mlat, nlev, & beglat, endlat, 1, nlev, ng_d, ng_d, 0) #else do k=1,nlev do j=beglat,endlat do i=1,nlon w3d(i,j,k) = q(i,j,k,n) enddo enddo enddo #endif if ( master ) then write(iuic) w3d call maxmin_global('Q', w3d, nlon, mlat, nlev, 1.) endif enddo if ( master .and. age_tracer ) write(iuic) age_time endif if( master ) then call close_file(iuic) write(*,*) 'FV restart file written and closed' endif !-------------------------------------------------------------- ! Write surface winds on A grid for coupler to restart properly !-------------------------------------------------------------- if ( master ) then !Modify the filename if it has an ensemble instance appendix call get_instance_filename('RESTART/fv_srf_wnd.res',restart_name ) restart_name=trim(restart_name) iuic = open_restart_file(restart_name, 'write') write(iuic) nlon, mlat endif #ifdef SPMD call mp_gather3d(ps_bp, w2d, nlon, mlat, 1, beglat, endlat, 1, 1, & 0, 0, 0) #endif if (master) write(iuic) w2d #ifdef SPMD call mp_gather3d(u_srf, w2d, nlon, mlat, 1, beglat, endlat, 1, 1, & 0, 0, 0) #endif if (master) write(iuic) w2d #ifdef SPMD call mp_gather3d(v_srf, w2d, nlon, mlat, 1, beglat, endlat, 1, 1, & 0, 0, 0) #endif if (master) then write(iuic) w2d call close_file(iuic) endif !---------------------------------------- ! Write additional tracers to netCDF file !---------------------------------------- if (do_fms_tracer_manager .and. ncnst > nt_phys) then !{ if(master) & & write(*,*)'Additional tracers present in field_table. These will be saved in a netcdf file.' #ifdef MARS_GCM call get_eta_level(nlev, p_ref, pfull, phalf, 0.01) #else call get_eta_level(nlev, 1.E5, pfull, phalf, 0.01) #endif MARS_GCM !Modify the filename if it has an ensemble instance appendix call get_instance_filename(fms_tracers_restart_file,restart_name ) restart_name=trim(restart_name) call mpp_open(iuic, restart_name, & action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE ) ! define axes call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', & data=lon ) call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', & data=lat ) call mpp_write_meta(iuic, axis_lev, 'pfull', 'mb', 'ref full pressure level', & data=pfull ) ! define fields allocate(tracer_fields(ncnst)) do n = nt_phys+1, ncnst call tr_get_tracer_names(MODEL_ATMOS, n, & tracer_name, tracer_longname, tracer_units) call mpp_write_meta(iuic, tracer_fields(n), & (/axis_lon, axis_lat, axis_lev/), & tracer_name, tracer_units, tracer_longname, pack=1) ! want doubles enddo ! write axes call mpp_write(iuic, axis_lon) call mpp_write(iuic, axis_lat) call mpp_write(iuic, axis_lev) ! write tracers do n = nt_phys+1, ncnst #ifdef SPMD call mp_gather3d(q(1,beglat-ng_d,1,n), w3d, nlon, mlat, nlev, & beglat, endlat, 1, nlev, ng_d, ng_d, 0) #else do k=1,nlev do j=beglat,endlat do i=1,nlon w3d(i,j,k) = q(i,j,k,n) enddo enddo enddo #endif if(master) call mpp_write(iuic, tracer_fields(n), data=w3d) enddo call mpp_close(iuic) deallocate(tracer_fields) endif !} deallocate ( w3d ) deallocate ( w2d ) else ! ! ------------------NetCDF restart files--------- call mpp_get_data_domain ( fv_domain, ids, ide, jds, jde ) call mpp_get_compute_domain( fv_domain, ics, ice, jcs, jce ) restart_name = trim(dflnm) ! add netcdf file suffix !Modify the filename if it has an ensemble instance appendix call get_instance_filename(restart_name,restart_name ) restart_name=trim(restart_name) call mpp_open(iuic, restart_name, & action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE ) ! define axes #ifdef MARS_GCM call get_eta_level(nlev, p_ref, pfull, phalf, 0.01) #else call get_eta_level(nlev, 1.E5, pfull, phalf, 0.01) #endif MARS_GCM ! if(master) write(*,*) "phalf values are ",phalf ! define longitude axis for u and tracer fields call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', & data=lon ) ! define longitude axis for v field call mpp_write_meta(iuic, axis_lonv, 'lonv', 'degree_E', 'longitude', & data=lonb(1:size(lonb)-1) ) ! define latitude axis for v and tracer fields call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', & data=lat ) ! define latitude axis for u field allocate(latu(size(lat))) ! The first latitude corresponds to meaningless data for the South Pole. ! For interpolation reasons, it has been decided to set this latitude to ! the second latitude with a small offset for Ferret compatibility. latu(1) = latb(2) - 1e-5 latu(2:) = latb(2:size(lat)) call mpp_write_meta(iuic, axis_latu, 'latu', 'degree_N', 'latitude', & data=latu ) ! define pressure axis call mpp_write_meta(iuic, axis_lev, 'pfull', 'mb', 'ref full pressure level', & data=pfull ) call mpp_write_meta(iuic, axis_levh, 'phalf', 'mb', 'ref half pressure level', & data=phalf ) call mpp_write_meta(iuic, axis_time, 'Time', 'time level', 'Time', & cartesian='T', data=dummy) call mpp_write_meta(iuic, ak_field, (/axis_levh/), & 'ak', 'Pa', 'ak-field', pack=1) ! want doubles call mpp_write_meta(iuic, bk_field, (/axis_levh/), & 'bk', ' ', 'bk-field', pack=1) call mpp_write_meta(iuic, 0, 'nlon', ival=nlon) call mpp_write_meta(iuic, 0, 'mlat', ival=mlat) call mpp_write_meta(iuic, 0, 'nlev', ival=nlev) call mpp_write_meta(iuic, 0, 'nt_phys', nt_phys) call mpp_write_meta(iuic, 0, 'days', days) call mpp_write_meta(iuic, 0, 'seconds', seconds) call mpp_write_meta(iuic, 0, 'ks', ks) call mpp_write_meta(iuic, u_field, (/axis_lon , axis_latu, axis_lev, axis_time/), & 'U', 'm/s', 'u-field', pack=1) call mpp_write_meta(iuic, v_field, (/axis_lonv, axis_lat , axis_lev, axis_time/), & 'V', 'm/s', 'v-field', pack=1) call mpp_write_meta(iuic, t_field, (/axis_lon , axis_lat , axis_lev, axis_time/), & 'T', 'K', 'Temperature', pack=1) call mpp_write_meta(iuic, phis_field, (/axis_lon, axis_lat, axis_time/), & 'Surface_geopotential', 'm**2/s**2', 'Surface Height', pack=1) call mpp_write_meta(iuic, delp_field, (/axis_lon, axis_lat, axis_lev, axis_time/), & 'DELP', 'pa', 'delp-field', pack=1) ! write axes call mpp_write(iuic, axis_time) call mpp_write(iuic, axis_lon) call mpp_write(iuic, axis_lonv) call mpp_write(iuic, axis_lat) call mpp_write(iuic, axis_latu) call mpp_write(iuic, axis_lev) call mpp_write(iuic, bk_field, bk) call mpp_write(iuic, ak_field, ak) call mpp_write(iuic, u_field, fv_domain, u(ids:ide,jds:jde,:)) call mpp_write(iuic, v_field, fv_domain, v(ids:ide,jds:jde,:)) call mpp_write(iuic, t_field, fv_domain, pt(ids:ide,jds:jde,:)) ! The mpp_io code uses Cray pointers which doesn't work well with the PSET ! 2D arrays of the FV code. Therefore a copy of the 2D array is needed before ! passing the array to the mpp_write routine. allocate (phis_alloc(ids:ide,jds:jde)) phis_alloc = -999.9 phis_alloc(ids:ide,jcs:jce) = phis(ids:ide,jcs:jce) call mpp_write(iuic, phis_field, fv_domain, phis_alloc) call mpp_write(iuic, delp_field, fv_domain, delp(ids:ide,jds:jde,:)) call mpp_close(iuic) call error_mesg ('write_fv_rst', 'Writing NetCDF formatted restart file: '//trim(restart_name), NOTE) call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) restart_name = trim(fms_tracers_restart_file)!//'1.nc' ! add netcdf file suffix !Modify the filename if it has an ensemble instance appendix call get_instance_filename(restart_name,restart_name ) restart_name=trim(restart_name) call mpp_open(iuic, restart_name, & action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE ) ! define axes call mpp_write_meta(iuic, axis_time, 'Time', 'time level', & 'Time',cartesian='T', data=dummy) ! define longitude axis for tracer fields call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', & data=lon ) ! define latitude axis for tracer fields call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', & data=lat ) ! define pressure axis call mpp_write_meta(iuic, axis_lev, 'pfull', 'mb', 'ref full pressure level', & data=pfull ) ! define fields allocate(tracer_fields(ntracers)) do n = 1,ntracers call tr_get_tracer_names(MODEL_ATMOS, n, & tracer_name, tracer_longname, tracer_units) call mpp_write_meta(iuic, tracer_fields(n), & (/axis_lon, axis_lat, axis_lev, axis_time/), & tracer_name, tracer_units, tracer_longname, pack=1) enddo ! write axes call mpp_write(iuic, axis_lon) call mpp_write(iuic, axis_lat) call mpp_write(iuic, axis_lev) call mpp_write(iuic, axis_time) do n = 1, ntracers call mpp_write(iuic, tracer_fields(n), fv_domain, q(ids:ide,jds:jde,:,n)) enddo deallocate(tracer_fields) call mpp_close(iuic) !Modify the filename if it has an ensemble instance appendix call get_instance_filename('RESTART/fv_srf_wnd.res',restart_name ) restart_name=trim(restart_name) call mpp_open(iuic, restart_name, & action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE ) ! define axes ! define longitude axis for u and tracer fields call mpp_write_meta(iuic, axis_lon, 'lon', 'degree_E', 'longitude', & data=lon ) ! define longitude axis for v field call mpp_write_meta(iuic, axis_lonv, 'lonv', 'degree_E', 'longitude', & data=lonb(1:size(lonb)-1) ) ! define latitude axis for v and tracer fields call mpp_write_meta(iuic, axis_lat, 'lat', 'degree_N', 'latitude', & data=lat ) ! define latitude axis for u field latu(1) = lat(1) latu(2:) = latb(2:size(lat)) call mpp_write_meta(iuic, axis_latu, 'latu', 'degree_N', 'latitude', & data=latu ) call mpp_write_meta(iuic, 'nlon', ival=nlon) call mpp_write_meta(iuic, 'mlat', ival=mlat) call mpp_write_meta(iuic, u_field, (/axis_lon , axis_latu/), & 'U_srf', 'm/s', 'surface u-field', pack=1) call mpp_write_meta(iuic, v_field, (/axis_lonv, axis_lat /), & 'V_srf', 'm/s', 'surface v-field', pack=1) call mpp_write_meta(iuic, phis_field, (/axis_lon, axis_lat/), & 'ps_bp', 'pa', 'surface pressure', pack=1) ! write axes call mpp_write(iuic, axis_lon) call mpp_write(iuic, axis_lonv) call mpp_write(iuic, axis_lat) call mpp_write(iuic, axis_latu) ! The mpp_io code uses Cray pointers which doesn't work well with the PSET ! 2D arrays of the FV code. Therefore a copy of the 2D array is needed before ! passing the array to the mpp_write routine. phis_alloc = -999.9 phis_alloc(ids:ide,jcs:jce) = ps_bp(ids:ide,jcs:jce) call mpp_write(iuic, phis_field, fv_domain, phis_alloc) phis_alloc(ids:ide,jcs:jce) = u_srf(ids:ide,jcs:jce) call mpp_write(iuic, u_field, fv_domain, phis_alloc) phis_alloc(ids:ide,jcs:jce) = v_srf(ids:ide,jcs:jce) call mpp_write(iuic, v_field, fv_domain, phis_alloc) endif ! end IF NETCDF call fv_array_sync() end subroutine write_fv_rst !--------------------------------------------------------------------- subroutine maxmin_global(qname, a, im, jm, km, fac) implicit none character*(*) qname integer im, jm, km real a(im,jm, km) real fac real pmax(km), pmin(km) real qmax, qmin integer i, j, k !$omp parallel do private(i, j, k) !Balaji: not parallelized, called with unshared arrays do k=1, km pmax(k) = a(1,1,k) pmin(k) = a(1,1,k) do j=1,jm do i=1, im pmax(k) = max (pmax(k), a(i,j,k)) pmin(k) = min (pmin(k), a(i,j,k)) enddo enddo enddo qmax = pmax(1) qmin = pmin(1) do k=2,km qmax = max (qmax, pmax(k)) qmin = min (qmin, pmin(k)) enddo write(6,*) qname, ' max = ', qmax * fac, ' min = ', qmin * fac end subroutine maxmin_global subroutine print_check_sum (label) #include "fv_arrays.h" character(len=*), intent(in) :: label integer :: n, chksum_unit integer :: isl, iel, jsl, jel integer(I8_KIND) :: zsum #ifdef use_shared_pointers #include "fv_point.inc" #endif call mpp_get_compute_domain( fv_domain, isl, iel, jsl, jel ) chksum_unit = stdout() if (mpp_pe() == mpp_root_pe()) write (chksum_unit,'(/a)') trim(label) ! sum all restart variables (except eta,peta) zsum = mpp_chksum(pt(isl:iel,jsl:jel,:)) if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(t)',zsum zsum = mpp_chksum(u(isl:iel,jsl:jel,:)) if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(u)',zsum zsum = mpp_chksum(v(isl:iel,jsl:jel,:)) if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(v)',zsum zsum = 0 ! sum of all tracers do n = 1, size(q,4) zsum = zsum + mpp_chksum(q(isl:iel,jsl:jel,:,n)) enddo if (mpp_pe() == mpp_root_pe()) write (chksum_unit,10) '(r)',zsum 10 format ('chksum',a7,' = ',z16) end subroutine print_check_sum end module fv_restart_mod