module my_esmf use esmf ! ESMF library module use params ,only: iulog use shr_kind_mod ,only: r8 => shr_kind_r8 use my_mpi ,only: ntask,ntaski,ntaskj,tasks,lon0,lon1,lat0,lat1,mytid,& nmagtaski,nmagtaskj,mlon0,mlon1,mlat0,mlat1 use apex ,only: gdlatdeg,gdlondeg use geogrid ,only: nlon,nlat,nlev,glon,glat,jspole,jnpole ! dynamically allocated geo grid use maggrid ,only: nmlev,nmlon,nmlonp1 use namelist ,only: model_name use timing ,only: time_esmf_init implicit none save public #include type(ESMF_Grid),save :: & geo_src_grid, mag_src_grid, & ! source grids (will not have periodic pts) geo_des_grid, mag_des_grid ! destination grids (will have periodic pts) ! ! 3d (i,j,k) ESMF Fields on geographic subdomains: ! type(ESMF_Field),save :: & ! 3d ESMF fields on geographic grid geo_ped, & ! pedersen conductivity geo_hal, & ! hall conductivity geo_zpot, & ! geopotential height (cm) geo_scht, & ! scale height (cm) geo_adotv1, & ! ue1 (m/s) geo_adotv2 ! ue2 (m/s) integer,parameter :: nf_3dgeo=6 ! number of 3d fields on geographic grid type(ESMF_Field),save :: f_3dgeo(nf_3dgeo) ! fields on 3d geo grid (could be bundled?) ! ! 2d (i,j) ESMF fields on geographic subdomains: ! type(ESMF_Field),save :: & ! 2d ESMF fields on geographic grid geo_sini, & ! sin(I_m) geo_adota1, & ! d(1)**2/D geo_adota2, & ! d(2)**2/D geo_a1dta2, & ! (d(1) dot d(2)) /D geo_be3 ! mag field strength (T) ! ! 3d (i,j,k) ESMF fields regridded to magnetic subdomains: ! type(ESMF_Field),save :: & ! 3d ESMF fields on geomagnetic grid mag_ped, & ! pedersen conductivity mag_hal, & ! hall conductivity mag_zpot, & ! geopotential height (cm) mag_scht, & ! scale height (cm) mag_adotv1, & ! ue1 (m/s) mag_adotv2 ! ue2 (m/s) ! ! 2d (i,j) ESMF fields on magnetic subdomains: ! type(ESMF_Field),save :: & ! 2d fields on geomagnetic grid mag_sini, & ! sin(I_m) mag_adota1, & ! d(1)**2/D mag_adota2, & ! d(2)**2/D mag_a1dta2, & ! (d(1) dot d(2)) /D mag_be3 ! mag field strength (T) ! ! 3d electric potential and electric field for mag to geo regridding: ! type(ESMF_Field),save :: mag_phi3d,mag_ephi3d,mag_elam3d,mag_emz3d type(ESMF_Field),save :: geo_phi3d,geo_ephi3d,geo_elam3d,geo_emz3d type(ESMF_RouteHandle),save :: & ! ESMF route handles for regridding routehandle_geo2mag, & ! for geo to mag regrid routehandle_mag2geo, & ! for mag to geo regrid routehandle_geo2mag_2d ! for 2d geo to mag ! ! Use either apex or esmf indices and weights for regridding: ! (in either case, ESMF SMM will be used for bilinear interpolation) ! Apex weights can be used only with TIMEGCM (not WACCM). If input ! model is WACCM, use_apex will be set false, and use_esmf will be true. ! logical :: use_apex = .true. logical :: use_esmf = .false. real(r8),save :: r8_nlon real(r8),allocatable,save :: unitv(:) ! private use_apex,use_esmf,routehandle_geo2mag, routehandle_mag2geo,& routehandle_geo2mag_2d,r8_nlon contains !----------------------------------------------------------------------- subroutine esmf_init ! ! This is called once per run. Initialize ESMF, create source and ! destination grids, set up weights for bilinear interpolation used ! in regridding. ! use apex,only: ig,jg,im,jm,wt,dim,djm ! ! Local: integer :: rc ! return code for ESMF calls integer :: k,j,i,n,ii,jj,indx,jndx,ijg real(ESMF_KIND_R8),pointer :: fptr(:,:,:) integer :: lbnd_destgeo(3),ubnd_destgeo(3) ! 3d bounds of destination geo grid integer :: lbnd_destmag(3),ubnd_destmag(3) ! 3d bounds of destination mag grid integer :: lbnd_srcgeo(3),ubnd_srcgeo(3) ! 3d bounds of source geo grid integer :: lbnd_srcmag(3),ubnd_srcmag(3) ! 3d bounds of source mag grid integer(ESMF_KIND_I4),pointer :: factorIndexList(:,:) ! esmf (2,nw) real(ESMF_KIND_R8),pointer :: factorList(:) ! esmf weights (nw) integer :: nw,nmaglats,nmaglons,ngeolats,ngeolons logical :: indxoob ! ! Apex indices and weights for geo2mag: integer(ESMF_KIND_I4),pointer :: igjg(:,:) ! apex geo indices for esmf (2,nw) real(ESMF_KIND_R8),pointer :: wts_geo2mag(:) ! apex weights for esmf (nw) ! ! Apex indices and weights for mag2geo: integer(ESMF_KIND_I4),pointer :: imjm(:,:) ! apex mag indices for esmf (2,nw) real(ESMF_KIND_R8),pointer :: wts_mag2geo(:) ! apex weights for esmf (nw) real(r8) :: starttime,endtime starttime = mpi_wtime() ! ! Initialize esmf: ! call ESMF_Initialize(logkindflag=ESMF_LOGKIND_MULTI,rc=rc) ! ! Set unit vector (this routine called once per run): ! allocate(unitv(nlon)) unitv(:) = 1._r8 r8_nlon = real(nlon) ! 8-byte float(nlon) ! ! Make magnetic and geographic grids for geo2mag regridding: ! call create_geo_grid(geo_src_grid,'src') ! geo source grid call create_mag_grid(mag_des_grid,'des') ! mag destination grid ! ! Make grids for mag2geo regridding: ! call create_mag_grid(mag_src_grid,'src') call create_geo_grid(geo_des_grid,'des') ! ! Create empty fields on geographic grid that will be transformed to ! the magnetic grid and passed as input to the dynamo. This does not ! assign any values. ! ! 3d fields on source geo grid (these exclude periodic points): ! call esmf_create_geofield(geo_ped,geo_src_grid, 'PED ',nlev) call esmf_create_geofield(geo_hal ,geo_src_grid, 'HAL ',nlev) call esmf_create_geofield(geo_zpot,geo_src_grid, 'ZPOT ',nlev) call esmf_create_geofield(geo_scht,geo_src_grid, 'SCHT ',nlev) call esmf_create_geofield(geo_adotv1,geo_src_grid,'ADOTV1 ',nlev) call esmf_create_geofield(geo_adotv2,geo_src_grid,'ADOTV2 ',nlev) ! ! Get 3d bounds of source geo field: ! call ESMF_FieldGet(geo_ped,localDe=0,farrayPtr=fptr, & computationalLBound=lbnd_srcgeo, & computationalUBound=ubnd_srcgeo,rc=rc) write(iulog,"('Bounds of source geo field: lbnd_srcgeo=',3i4,' ubnd_srcgeo=',3i4)") & lbnd_srcgeo,ubnd_srcgeo ! ! 2d fields on source geo grid (these exclude periodic points): ! call esmf_create_geofield(geo_sini ,geo_src_grid,'SINI ',0) call esmf_create_geofield(geo_adota1,geo_src_grid,'ADOTA1 ',0) call esmf_create_geofield(geo_adota2,geo_src_grid,'ADOTA2 ',0) call esmf_create_geofield(geo_a1dta2,geo_src_grid,'A1DTA2 ',0) call esmf_create_geofield(geo_be3 ,geo_src_grid,'BE3 ',0) ! ! 3d fields on destination mag grid (will include periodic point): ! call esmf_create_magfield(mag_ped ,mag_des_grid, 'PED ',nmlev) call esmf_create_magfield(mag_hal ,mag_des_grid, 'HAL ',nmlev) call esmf_create_magfield(mag_zpot,mag_des_grid, 'ZPOT ',nmlev) call esmf_create_magfield(mag_scht,mag_des_grid, 'SCHT ',nmlev) call esmf_create_magfield(mag_adotv1,mag_des_grid,'ADOTV1 ',nmlev) call esmf_create_magfield(mag_adotv2,mag_des_grid,'ADOTV2 ',nmlev) ! ! Get 3d bounds of destination mag field: ! call ESMF_FieldGet(mag_ped,localDe=0,farrayPtr=fptr, & computationalLBound=lbnd_destmag, & computationalUBound=ubnd_destmag,rc=rc) write(iulog,"('Bounds of destination mag field: lbnd_destmag=',3i4,' ubnd_destmag=',3i4)") & lbnd_destmag,ubnd_destmag ! ! 2d fields on destination mag grid (will include periodic point): ! call esmf_create_magfield(mag_sini ,mag_des_grid,'SINI ',0) call esmf_create_magfield(mag_adota1,mag_des_grid,'ADOTA1 ',0) call esmf_create_magfield(mag_adota2,mag_des_grid,'ADOTA2 ',0) call esmf_create_magfield(mag_a1dta2,mag_des_grid,'A1DTA2 ',0) call esmf_create_magfield(mag_be3 ,mag_des_grid,'BE3 ',0) ! ! 3d fields on source mag grid for mag2geo: ! call esmf_create_magfield(mag_phi3d ,mag_src_grid,'PHIM3D ',nmlev) call esmf_create_magfield(mag_ephi3d,mag_src_grid,'EPHI3D ',nmlev) call esmf_create_magfield(mag_elam3d,mag_src_grid,'ELAM3D ',nmlev) call esmf_create_magfield(mag_emz3d ,mag_src_grid,'EMZ3D ',nmlev) ! ! 3d fields on destination geo grid for mag2geo: ! call esmf_create_geofield(geo_phi3d ,geo_des_grid,'PHIG3D ',nlev) call esmf_create_geofield(geo_ephi3d,geo_des_grid,'EPHI3D ',nlev) call esmf_create_geofield(geo_elam3d,geo_des_grid,'ELAM3D ',nlev) call esmf_create_geofield(geo_emz3d ,geo_des_grid,'EMZ3D ',nlev) ! ! Get 3d bounds of destination geo field: ! call ESMF_FieldGet(geo_phi3d,localDe=0,farrayPtr=fptr,& computationalLBound=lbnd_destgeo, & computationalUBound=ubnd_destgeo,rc=rc) ! ! Use either apex or esmf weights in esmf SMM regridding: ! (only esmf weights are allowed with WACCM input model) ! if (trim(model_name)=='WACCM'.and.use_apex) then write(iulog,"('>>> WARNING: cannot use apex weights with ',a,' input model.')") & model_name write(iulog,"('Will use esmf indices and weights')") use_apex = .false. use_esmf = .true. endif if ((use_apex.and.use_esmf).or. & (.not.use_apex.and..not.use_esmf)) then write(iulog,"(a,a,l1,a,l1)") & '>>> Must choose between use_apex and use_esmf:',& ' use_apex=',use_apex,' use_esmf=',use_esmf call shutdown('use_apex or use_esmf') endif ! ! Save route handles for grid transformations in both directions ! geo2mag and mag2geo. FieldRegridStore needs to be called only ! once for each transformation before the timestep loop (src and ! dest fields are still required, so just use ped here). Once inside ! the timestep loop, the same routehandle can be used for all fields ! that are regridded in the given direction. ! ! These calls will leave *.vtk info files in execdir: ! call ESMF_GridWriteVTK(geo_src_grid, & ! staggerloc=ESMF_STAGGERLOC_CENTER, filename="geoGrid",rc=rc) ! call ESMF_GridWriteVTK(mag_des_grid, & ! staggerloc=ESMF_STAGGERLOC_CENTER, filename="magGrid",rc=rc) ! ! Save route handle and get esmf indices and weights for geo2mag: ! ! 0: >>> esmf_init: error return from ESMF_FieldRegridStore for 3d geo2mag: rc= 513 ! PET0 ESMF_FieldRegrid.F90:771 ESMF_FieldRegridStoreNX Arguments are incompatible ! call ESMF_FieldRegridStore(srcField=geo_ped,dstField=mag_ped, & regridMethod=ESMF_REGRIDMETHOD_BILINEAR, & routeHandle=routehandle_geo2mag,factorIndexList=factorIndexList, & factorList=factorList,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(a,a,i4)") '>>> esmf_init: error return from ', & 'ESMF_FieldRegridStore for 3d geo2mag: rc=',rc call shutdown('ESMF_FieldRegridStore') endif ! ! Set up apex indices and weights for geo2mag (TIMEGCM only): ! if (use_apex) then ! Not with WACCM input file ! ! number of mag lats in current subdomain (mlat1-mlat0+1) ! number of mag lons in current subdomain (mlon1-mlon0+1) ! nmaglats = tasks(mytid)%nmaglats nmaglons = tasks(mytid)%nmaglons ! ! nw = number of weights = 4*nmaglons*nmaglats ! This is all 4 corners of box at each mag subdomain point i,j ! size(factorList) should be same as 4*nmaglons*nmaglats ! nw = 4*nmaglons*nmaglats if (size(factorList) /= 4*nmaglons*nmaglats) then write(iulog,"('>>> size(factorList)=',i5,' /= 4*nmaglons*nmaglats=',i5)") & size(factorList),4*nmaglons*nmaglats call shutdown('size of factorList for mag2geo') endif ! ! Apex weights and indices on magnetic subdomain, to be passed to esmf ! for SMM regridding: ! allocate(wts_geo2mag(nw)) ! weights (factorList) allocate(igjg(2,nw)) ! indices (factorIndexList) ! ! Loop over destination mag subdomain and n-box: k = 1 do j=1,nmaglats do i=1,nmaglons do n=1,4 ! ! Global indices into local subdomain: ! (mlon0,1 and mlat0,1 match lbnd,ubnd of dest mag grid subdomain) ii = mlon0+(i-1) ! i index into ig,jg jj = mlat0+(j-1) ! j index into ig,jg igjg(2,k) = (jj-1)*nmlonp1+ii ! seq index for mag dest field ! ! indx = i index to nth corner, jndx = j index to nth corner ! Init to lower-left corner (ig,jg): ! indx = ig(ii,jj) jndx = jg(ii,jj) if (n==2.or.n==4) indx = indx+1 if (n==3.or.n==4) jndx = jndx+1 indxoob = .false. if (indx > nlon) then ! out of bounds, so wrap to i=1 indx = 1 indxoob = .true. endif ! ! ijg = sequence index to nth corner of src geo grid: ijg = indx+(jndx-1)*nlon ! global seq index into geo src grid igjg(1,k) = ijg ! apex factor index list for src grid wts_geo2mag(k) = wt(n,ii,jj) ! apex weight at n ! ! apex uses a counterclockwise increment of n, whereas esmf uses left-to-right. ! Swap n = 3 and 4 to accomodate apex ordering: if (n == 3) wts_geo2mag(k) = wt(4,ii,jj) ! swap 3 and 4 for apex order if (n == 4) wts_geo2mag(k) = wt(3,ii,jj) if (indxoob) then igjg(1,k) = igjg(1,k-1) igjg(1,k-1) = ijg wts_geo2mag(k) = wts_geo2mag(k-1) wts_geo2mag(k-1) = wt(n,ii,jj) if (n == 3) wts_geo2mag(k-1) = wt(4,ii,jj) ! swap 3 and 4 for apex order if (n == 4) wts_geo2mag(k-1) = wt(3,ii,jj) endif ! ! Increment sequence index in igjg(1:2,k) and wts_geo2mag(k) k = k+1 enddo ! do n=1,4 enddo ! do i=1,nmaglons enddo ! do j=1,nmaglats ! ! End if using apex indices and weights for mag2geo endif ! use_apex ! ! Store route handle for geo2mag 3d fields: ! ! Use apex index list and weights for 3d fields: ! 2: >>> esmf_init: error return from ESMF_FieldSMMStore for 3d geo2mag: rc= 506 ! 3: >>> esmf_init: error return from ESMF_FieldSMMStore for 3d geo2mag: rc= 506 ! PET2 ESMCI_Array.C:7214 ESMCI::Array::sparseMatMulStore ! ( Invalid argument - factorIndexList contains srcSeqIndex outside srcArray bounds ! if (use_apex) then write(iulog,"('Using apex indices and weights in FieldSMMStore for 3d geo2mag')") call ESMF_FieldSMMStore(geo_ped,mag_ped,routehandle_geo2mag, & wts_geo2mag,igjg,rc=rc) else ! ! OR use esmf index list and weights for 3d fields: ! write(iulog,"('Using esmf indices and weights in FieldSMMStore for 3d geo2mag')") call ESMF_FieldSMMStore(geo_ped,mag_ped,routehandle_geo2mag, & factorList,factorIndexList,rc=rc) endif if (rc /= ESMF_SUCCESS) then write(iulog,"(2a,i4)") '>>> esmf_init: error return from ESMF_FieldSMMStore for ',& '3d geo2mag: rc=',rc call shutdown('ESMF_FieldSMMStore') endif ! ! Store route handle geo2mag 2d fields: ! ! Use apex index list and weights for 2d fields: if (use_apex) then write(iulog,"('Using apex indices and weights in FieldSMMStore for 2d geo2mag')") call ESMF_FieldSMMStore(geo_sini,mag_sini,routehandle_geo2mag_2d, & wts_geo2mag,igjg,rc=rc) else ! ! OR use esmf index list and weights for 2d fields: write(iulog,"('Using esmf indices and weights in FieldSMMStore for 2d geo2mag')") call ESMF_FieldSMMStore(geo_sini,mag_sini,routehandle_geo2mag_2d, & factorList,factorIndexList,rc=rc) endif ! if (rc /= ESMF_SUCCESS) then write(iulog,"(2a,i4)") '>>> esmf_init: error return from ESMF_FieldSMMStore',& ' for 2d geo2mag: rc=',rc call shutdown('esmf_FieldRegridStore') endif ! ! Save route handle and get esmf indices and weights for mag2geo: ! (this overwrites factorIndexList and factorList from geo2mag call above) ! call ESMF_FieldRegridStore(srcField=mag_phi3d,dstField=geo_phi3d, & regridMethod=ESMF_REGRIDMETHOD_BILINEAR, & routeHandle=routehandle_mag2geo,factorIndexList=factorIndexList,& factorList=factorList,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(2a,i4)") '>>> esmf_init: error return from ',& 'ESMF_FieldRegridStore for 3d mag2geo: rc=',rc call shutdown('esmf_FieldRegridStore') endif ! ! Set up apex weights and indices for mag2geo: ! (WACCM has geographic poles, TIMEGCM does not, but cannot use apex ! weights with WACCM, so it doesnt really matter) ! if (use_apex) then ! Not with WACCM input file ngeolats = tasks(mytid)%nlats do j=lat0,lat1 if ((j==1.and.glat(j)/=-90._r8).or.(j==nlat.and.glat(j)/=90._r8)) & ngeolats = ngeolats+1 ! 0 and 37 for poles (TIMEGCM only) enddo ngeolons = tasks(mytid)%nlons nw = 4*ngeolons*ngeolats if (size(factorList) /= 4*ngeolons*ngeolats) then write(iulog,"('>>> size(factorList)=',i5,' /= 4*ngeolons*ngeolats=',i5)") & size(factorList),4*ngeolons*ngeolats write(iulog,"('tasks(mytid)%nlats=',i4,' ngeolats=',i4,' ngeolons=',i4)") & tasks(mytid)%nlats,ngeolats,ngeolons call shutdown('wrong nw for mag2geo') endif if (ubnd_destgeo(2)-lbnd_destgeo(2)+1 /= ngeolats) then write(iulog,"(2a,i4,a,i4,a,i4)") '>>> dest geo bounds /= ngeolats: ', & 'ubnd_destgeo(2)=',ubnd_destgeo(2),' lbnd_destgeo(2)=',lbnd_destgeo(2),& ' ngeolats=',ngeolats call shutdown('dest geo bounds /= ngeolats') endif if (ubnd_destgeo(1)-lbnd_destgeo(1)+1 /= ngeolons) then write(iulog,"(2a,i4,a,i4,a,i4)") '>>> dest geo bounds /= ngeolons: ', & 'ubnd_destgeo(1)=',ubnd_destgeo(1),' lbnd_destgeo(1)=',lbnd_destgeo(1),& ' ngeolats=',ngeolats call shutdown('dest geo bounds /= ngeolons') endif allocate(wts_mag2geo(nw)) ! weights (factorList) allocate(imjm(2,nw)) ! indices (factorIndexList) for src,dest grids ! ! Loop over destination geo grid subdomain and n-box, using upper and ! lower bounds of destination geo grid obtained from ESMF_FieldGet call above. ! For example, if ntask=4, j=0,18 or j=19,37, i= 1,38 or i=39,76 ! k = 1 do j=lbnd_destgeo(2),ubnd_destgeo(2) ! latitude bounds do i=lbnd_destgeo(1),ubnd_destgeo(1) ! longitude bounds do n=1,4 ! global seq index for geo dest field (WARNING: nlon+4 is for TIMEGCM only) imjm(2,k) = j*(nlon+4)+i ! ! ii,jj are geo global indices into apex im,jm and dim,djm: ! Sort out longitude end points (note longitude dimension is nlonp1) ! ii = i-2 ! i=3,74 -> ii=1,72 if (i <= 2) ii = nlon-2+i ! i=1,2 -> ii=71,72 if (i >= nlon+3) ii = i-nlon-2 ! i=75,76 -> ii=1,2 jj = j indx = im(ii,jj) jndx = jm(ii,jj) if (n==2.or.n==4) indx = indx+1 if (n==3.or.n==4) jndx = jndx+1 ! ! Check for magnetic longitude end points: ! When indx==80, then n is always 2 or 4 ! When indx== 1, then n is always 1 or 3 ! indxoob = .false. if (indx == nmlon.and.(n==1.or.n==3)) then indx = 1 indxoob = .true. endif if (indx > nmlon.and.(n==2.or.n==4)) then indx = nmlon indxoob = .true. endif ! ! ijg = sequence index to nth corner of src geo grid: ijg = indx+(jndx-1)*nmlon ! global seq index into mag src grid imjm(1,k) = ijg ! apex factor index list for mag grid !my_esmf.F90(488): error #6633: The type of the actual argument differs from the type of the dummy argument. [DIM] wts_mag2geo(k) = select_wt_mag2geo(n,dim(ii,jj),djm(ii,jj)) ! ! Sort out consistent ordering of corner weights (apex ~ esmf): if (indxoob) then if (n==1) wts_mag2geo(k) = & select_wt_mag2geo(2,dim(ii,jj),djm(ii,jj)) if (n==2) wts_mag2geo(k) = & select_wt_mag2geo(1,dim(ii,jj),djm(ii,jj)) else if (n == 3) wts_mag2geo(k) = & select_wt_mag2geo(4,dim(ii,jj),djm(ii,jj)) if (n == 4) wts_mag2geo(k) = & select_wt_mag2geo(3,dim(ii,jj),djm(ii,jj)) endif k = k+1 enddo ! n=1,4 enddo ! i-loop over geo subdomain longitudes enddo ! j-loop over geo subdomain latitudes ! ! End if using apex indices and weights for geo2mag: endif ! use_apex ! ! mag2geo 3d fields: ! ! Use apex index list and weights for 3d fields: if (use_apex) then write(iulog,"('Using apex indices and weights in FieldSMMStore for 3d mag2geo')") call ESMF_FieldSMMStore(mag_phi3d,geo_phi3d,routehandle_mag2geo,& wts_mag2geo,imjm,rc=rc) else ! ! OR use esmf index list and weights for 3d fields: ! write(iulog,"('Using esmf indices and weights in FieldSMMStore for 3d mag2geo')") call ESMF_FieldSMMStore(mag_phi3d,geo_phi3d,routehandle_mag2geo,& factorList,factorIndexList,rc=rc) endif if (rc /= ESMF_SUCCESS) then write(iulog,"(2a,i4)") '>>> esmf_init: error return from ESMF_FieldSMMStore ',& 'for 3d mag2geo: rc=',rc call shutdown('ESMF_FieldSMMStore') endif endtime = mpi_wtime() time_esmf_init=time_esmf_init+(endtime-starttime) end subroutine esmf_init !----------------------------------------------------------------------- real(r8) function select_wt_mag2geo(n,dimx,djmy) integer,intent(in) :: n real(r8),intent(in) :: dimx,djmy select_wt_mag2geo = 0._r8 select case (n) case(1) select_wt_mag2geo = (1._r8-dimx)*(1._r8-djmy) case(2) select_wt_mag2geo = dimx*(1._r8-djmy) case(3) select_wt_mag2geo = dimx*djmy case(4) select_wt_mag2geo = (1._r8-dimx)*djmy end select end function select_wt_mag2geo !----------------------------------------------------------------------- subroutine create_geo_grid(grid_out,srcdes) ! ! Args: type(ESMF_Grid),intent(out) :: grid_out character(len=*),intent(in) :: srcdes ! ! Local: integer :: i,j,n,rc integer :: lbnd_lat,ubnd_lat,lbnd_lon,ubnd_lon,lbnd(1),ubnd(1) real(ESMF_KIND_R8),pointer :: coordX(:),coordY(:) integer :: nlons_task(ntaski) ! number of lons per task integer :: nlats_task(ntaskj) ! number of lats per task logical :: has_poles ! ! We are creating either a source grid or a destination grid: ! if (srcdes /= 'src' .and. srcdes /= 'des') then write(iulog,"(a)") '>>> create_mag_grid: srcdes = ''',srcdes, & ''' but must be either ''src'' or ''des''' call shutdown('srcdes') endif ! ! nlons_task(ntaski) = number of geo lons per task. ! do i=1,ntaski loop: do n=0,ntask-1 if (tasks(n)%mytidi==i-1) then nlons_task(i) = tasks(n)%nlons exit loop endif enddo loop enddo ! ! Exclude periodic points (2 points fewer for procs at each end) ! for source grids only (east and west edges of task table). ! (TIMEGCM only) ! if (srcdes == 'src'.and.trim(model_name)=='TIMEGCM') then do n=0,ntask-1 if (tasks(n)%mytidi==ntaski-1.or.tasks(n)%mytidi==0) & nlons_task(tasks(n)%mytidi+1) = tasks(n)%nlons-2 enddo endif ! ! nlats_task(ntaskj) = number of geo lats per task. ! do j=1,ntaskj loop1: do n=0,ntask-1 if (tasks(n)%mytidj==j-1) then nlats_task(j) = tasks(n)%nlats exit loop1 endif enddo loop1 enddo ! ! Check to see if global glat(nlat) has poles (WACCM does, TIMEGCM does not): has_poles = .false. do j=1,nlat if (abs(glat(j))==90._r8) has_poles = .true. enddo ! write(iulog,"('create_geo_grid: model=',a,' has_poles=',l1)") model_name,has_poles ! ! If glat does not have poles, add extra points at north and south edges ! of task table (source or destination grid): ! if (.not.has_poles) then ! probably TIMEGCM do n=0,ntask-1 if (tasks(n)%mytidj==ntaskj-1.or.tasks(n)%mytidj==0) & nlats_task(tasks(n)%mytidj+1) = tasks(n)%nlats+1 enddo ! ! Create 2d geographic grid (minimum lat index is 0 to include poles): grid_out = ESMF_GridCreate1PeriDim( & countsPerDEDim1=nlons_task, coordDep1=(/1/), & countsPerDEDim2=nlats_task, coordDep2=(/2/), & indexflag=ESMF_INDEX_GLOBAL,minIndex=(/1,0/),rc=rc) else ! glat does have poles ! ! Create 2d geographic grid (grid already has poles): grid_out = ESMF_GridCreate1PeriDim( & countsPerDEDim1=nlons_task, coordDep1=(/1/), & countsPerDEDim2=nlats_task, coordDep2=(/2/), & indexflag=ESMF_INDEX_GLOBAL,minIndex=(/1,1/),rc=rc) endif if (rc /=ESMF_SUCCESS) then write(iulog,"(/,2a,i4)") '>>> create_geo_grid: error return from ',& 'ESMF_GridCreate1PeriDim: rc=',rc call shutdown('create_geo_grid') endif ! ! Allocate coordinates: ! call ESMF_GridAddCoord(grid_out,staggerloc=ESMF_STAGGERLOC_CENTER,rc=rc) if (rc /=ESMF_SUCCESS) then write(iulog,"(/,a)") '>>> create_geo_grid: error return from ESMF_GridAddCoord' call shutdown('create_geo_grid') endif ! ! Get pointer and set geo grid longitude coordinates: ! call ESMF_GridGetCoord(grid_out, coordDim=1, localDE=0, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordX, staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /=ESMF_SUCCESS) then write(iulog,"(/,2a)") '>>> create_geo_grid: error return from ',& 'ESMF_GridGetCoord for longitude coords' call shutdown('create_geo_grid') endif lbnd_lon = lbnd(1) ; ubnd_lon = ubnd(1) if (srcdes == 'src') then ! source grid excludes periodic points do i=lbnd_lon,ubnd_lon coordX(i) = glon(i) ! 1 -> 72 enddo else ! destination grid includes periodic points do i=lbnd_lon,ubnd_lon if (i <= 2) then ! comments assume dlon = 5 deg coordX(i) = glon(nlon-2+i) ! i=1,2 <- 73,74 == glon(nlon-1,nlon) elseif (i >= nlon+3) then coordX(i) = glon(i-nlon-2) ! i=75,76 <- 3,4 == glon(1,2) else coordX(i) = glon(i-2) ! i=3,74 == glon(1,nlon) endif enddo endif ! ! Get pointer and set geo grid latitude coordinates, including poles: ! call ESMF_GridGetCoord(grid_out, coordDim=2, localDE=0, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordY, staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /=ESMF_SUCCESS) then write(iulog,"(/,2a)") '>>> create_geo_grid: error return from ',& 'ESMF_GridGetCoord for latitude coords' call shutdown('create_geo_grid') endif lbnd_lat = lbnd(1) ; ubnd_lat = ubnd(1) do j=lbnd_lat,ubnd_lat if (j==jspole) then coordY(j) = -90._r8 elseif (j==jnpole) then coordY(j) = +90._r8 else coordY(j) = glat(j) endif enddo write(iulog,"(4a,2i4,a,2i4,a,2i4,a,2i4)") 'Created ESMF ',srcdes,' geo grid: ', & ' lbnd,ubnd_lon=',lbnd_lon,ubnd_lon,' lon0,1=',lon0,lon1,& ' lbnd,ubnd_lat=',lbnd_lat,ubnd_lat,' lat0,1=',lat0,lat1 end subroutine create_geo_grid !----------------------------------------------------------------------- subroutine esmf_create_geofield(field,grid,name,nlev) ! ! Create ESMF field (2d or 3d) on geo grid (will exclude periodic points) ! If nlev == 0, field is 2d (i,j), otherwise field is 3d, ! and 3rd dimension is ungridded ! ! Args: integer,intent(in) :: nlev ! if nlev == 0, field is 2d (i,j) type(ESMF_Grid),intent(in) :: grid character(len=*),intent(in) :: name type(ESMF_Field),intent(out) :: field ! ! Local: integer :: rc type(ESMF_ArraySpec) :: arrayspec ! ! Create 3d field (i,j,k), with non-distributed vertical dimension: if (nlev > 0) then call ESMF_ArraySpecSet(arrayspec,3,ESMF_TYPEKIND_R8,rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_ArraySpecSet') field = ESMF_FieldCreate(grid, arrayspec,ungriddedLBound=(/1/), & ungriddedUBound=(/nlev/),staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_FieldCreate 3d') ! ! Create 2d field (i,j): else ! create 2d field call ESMF_ArraySpecSet(arrayspec,2,ESMF_TYPEKIND_R8,rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_ArraySpecSet') field = ESMF_FieldCreate(grid, arrayspec,& staggerloc=ESMF_STAGGERLOC_CENTER,rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_FieldCreate 2d') endif end subroutine esmf_create_geofield !----------------------------------------------------------------------- subroutine create_mag_grid(grid_out,srcdes) ! ! Args: type(ESMF_Grid),intent(out) :: grid_out character(len=*),intent(in) :: srcdes ! ! Local: integer :: istat,i,j,n,rc real(ESMF_KIND_R8),pointer :: coordX(:,:),coordY(:,:) integer :: lbnd(2),ubnd(2) integer :: nmlons_task(ntaski) ! number of lons per task integer :: nmlats_task(ntaskj) ! number of lats per task ! ! We are creating either a source grid or a destination grid: ! if (srcdes /= 'src' .and. srcdes /= 'des') then write(iulog,"(a)") '>>> create_mag_grid: srcdes = ''',srcdes, & ''' but must be either ''src'' or ''des''' call shutdown('srcdes') endif ! ! nmlons_task(nmagtaski) = number of mag lons per task in lon dim ! do i=1,nmagtaski loop: do n=0,ntask-1 if (tasks(n)%magtidi==i-1) then nmlons_task(i) = tasks(n)%nmaglons exit loop endif enddo loop enddo ! ! Exclude periodic points (1 point fewer for mpi tasks at east end) ! for source grids (this overwrites above for eastern-most tasks): ! if (srcdes == 'src') then do n=0,ntask-1 if (tasks(n)%magtidi==nmagtaski-1) then ! east edge of proc matrix nmlons_task(tasks(n)%magtidi+1) = tasks(n)%nmaglons-1 endif enddo endif ! ! nmlats_task(nmagtaskj) = number of mag lats per task in lat dim ! do j=1,nmagtaskj loop1: do n=0,ntask-1 if (tasks(n)%magtidj==j-1) then nmlats_task(j) = tasks(n)%nmaglats exit loop1 endif enddo loop1 enddo ! ! Create curvilinear magnetic grid (both coords depend ! on both dimensions, i.e., lon(i,j),lat(i,j)): ! grid_out = ESMF_GridCreate1PeriDim( & countsPerDEDim1=nmlons_task, coordDep1=(/1,2/), & countsPerDEDim2=nmlats_task, coordDep2=(/1,2/), & indexflag=ESMF_INDEX_GLOBAL,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(2a,i4)") '>>> create_mag_grid: error return from ',& 'ESMF_GridCreateShapeTile: rc=',rc call shutdown('create_mag_grid') endif ! ! Allocate coordinates: ! call ESMF_GridAddCoord(grid_out,staggerloc=ESMF_STAGGERLOC_CENTER,rc=rc) if (rc /=ESMF_SUCCESS) then write(iulog,"(2a,i4)") '>>> create_mag_grid: error return from ',& 'ESMF_GridAddCoord: rc=',rc call shutdown('create_mag_grid') endif ! ! Get pointer and set mag grid longitude coordinates: ! call ESMF_GridGetCoord(grid_out, coordDim=1, localDE=0, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordX, staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(i4)") '>>> create_mag_grid: error return from ', & 'ESMF_GridGetCoord for longitude coords: rc=',rc call shutdown('create_mag_grid') endif do j=lbnd(2),ubnd(2) do i=lbnd(1),ubnd(1) coordX(i,j) = gdlondeg(i,j) enddo enddo ! ! Get pointer and set mag grid latitude coordinates: ! call ESMF_GridGetCoord(grid_out, coordDim=2, localDE=0, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordY, staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(i4)") '>>> create_mag_grid: error return from ',& 'ESMF_GridGetCoord for latitude coords: rc=',rc call shutdown('create_mag_grid') endif do j=lbnd(2),ubnd(2) do i=lbnd(1),ubnd(1) coordY(i,j) = gdlatdeg(i,j) enddo enddo write(iulog,"(4a,2i4,a,2i4,a,2i4,a,2i4)") 'Created ESMF ',srcdes,' mag grid: ', & ' lbnd,ubnd_lon=',lbnd(1),ubnd(1),' mlon0,1=',mlon0,mlon1, & ' lbnd,ubnd_lat=',lbnd(2),ubnd(2),' mlat0,1=',mlat0,mlat1 end subroutine create_mag_grid !----------------------------------------------------------------------- subroutine esmf_create_magfield(field,grid,name,nlev) ! ! Create ESMF field (2d or 3d) on mag grid. This will include the ! mag periodic point, which will be zero after regridding. ! If nlev == 0, field is 2d (i,j), otherwise field is 3d, ! and 3rd dimension is ungridded ! ! Args: integer,intent(in) :: nlev ! if nlev == 0, field is 2d (i,j) type(ESMF_Grid),intent(in) :: grid character(len=*),intent(in) :: name type(ESMF_Field),intent(out) :: field ! ! Local: integer :: rc type(ESMF_ArraySpec) :: arrayspec type(ESMF_Array) :: array3d,array2d type(ESMF_DistGrid) :: distgrid ! ! Get necessary information from the mag grid: call ESMF_GridGet(grid, staggerloc=ESMF_STAGGERLOC_CENTER,& distgrid=distgrid,rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_GridGet') ! ! Create 3d mag field (i,j,k), with non-distributed vertical dimension: ! (add periodic point in longitude with computationalEdgeUWidth) ! if (nlev > 0) then call ESMF_ArraySpecSet(arrayspec,3,ESMF_TYPEKIND_R8,rc=rc) if (rc /= ESMF_SUCCESS)call shutdown('ESMF_ArraySpecSet 3d mag') array3d = ESMF_ArrayCreate(arrayspec=arrayspec, & distgrid=distgrid,computationalEdgeUWidth=(/1,0/), & undistLBound=(/1/),undistUBound=(/nlev/), & indexflag=ESMF_INDEX_GLOBAL,rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_ArrayCreate 3d mag') field = ESMF_FieldCreate(grid, array3d, & ungriddedLBound=(/1/), ungriddedUBound=(/nlev/), & staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_FieldCreate 3d mag') ! ! Create 2d mag field (i,j): ! (add periodic point in longitude with computationalEdgeUWidth) ! else ! create 2d field call ESMF_ArraySpecSet(arrayspec,2,ESMF_TYPEKIND_R8,rc=rc) if (rc /= ESMF_SUCCESS)call shutdown('ESMF_ArraySpecSet 2d mag') array2d = ESMF_ArrayCreate(arrayspec=arrayspec, & distgrid=distgrid,computationalEdgeUWidth=(/1,0/), & indexflag=ESMF_INDEX_GLOBAL,rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_ArrayCreate 2d mag') field = ESMF_FieldCreate(grid, array2d, & staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (rc /= ESMF_SUCCESS) call shutdown('ESMF_FieldCreate 2d mag') endif end subroutine esmf_create_magfield !----------------------------------------------------------------------- subroutine esmf_set2d_geo(field,grid,fname,f,ilon0,ilon1,ilat0,ilat1) use my_mpi,only: mp_geopole_2d ! ! Set values of a 2d ESMF field on geographic source grid, prior to ! geographic to magnetic grid transformation. (Essentially the same ! as esmf_set3d_geo, except for 2d fields instead of 3d) ! Periodic points are excluded, geographic poles are at j==jspole and jnpole ! ! Args: type(ESMF_Field) ,intent(in) :: field type(ESMF_Grid) ,intent(in) :: grid character(len=*) ,intent(in) :: fname ! field name integer ,intent(in) :: ilon0,ilon1,ilat0,ilat1 real(r8) ,intent(in) :: f(ilon0:ilon1,ilat0:ilat1) ! ! Local: integer :: i,ii,j,rc real(ESMF_KIND_R8),pointer :: fptr(:,:) ! i,j integer :: lbnd(2),ubnd(2) ! ! fpole_jpm2 is for call to mp_geopole_2d for 1 field only. When this ! is working, mp_geopole_2d should be called for all 2d fields before ! this routine is called. ! real(r8) :: fpole_jpm2(nlon,4,1) ! 1 field only ! ! Get pointer to the field: call ESMF_FieldGet(field,localDe=0,farrayPtr=fptr,& computationalLBound=lbnd,computationalUBound=ubnd,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(a,i4)") '>>> esmf_set2d_geo: error from ESMF_FieldGet: rc=',rc call shutdown('esmf_set2d_geo:ESMF_FieldGet') endif ! fptr(:,:) = 0._r8 ! init ! ! Set interior latitudes (excluding poles): do j=lbnd(2),ubnd(2) if (j /= jspole .and. j /= jnpole) then do i=lbnd(1),ubnd(1) ii = i if (trim(model_name)=='TIMEGCM') & ii = i+2 ! timegcm data is at 3->74 fptr(i,j) = f(ii,j) enddo endif ! interior latitudes only enddo ! ! Call mp_geopole_2d to get global lons surrounding the 2 lats ! adjacent to each pole. ! subroutine mp_geopole_2d(f,ilon0,ilon1,ilat0,ilat1, ! nglblon,jspole,jnpole,fpole_jpm2,nf) ! Return fpole_jpm2(nglblon,1->4,nf) as: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) ! call mp_geopole_2d(fptr,lbnd(1),ubnd(1),lbnd(2),ubnd(2),& nlon,jspole,jnpole,fpole_jpm2,1) ! ! Set fptr poles using global lons from mp_geopole_2d: do j=lbnd(2),ubnd(2) ! lat ! ! Set south pole: if (j==jspole) then ! south pole if (lbnd(2) > jspole+2) & ! need lats 1,2 for south pole call shutdown('esmf_set2d_geo need j=1,2') fptr(lbnd(1),j) = & (9._r8*dot_product(unitv,fpole_jpm2(:,1,1))- & dot_product(unitv,fpole_jpm2(:,2,1))) & / (8._r8*r8_nlon) do i=lbnd(1)+1,ubnd(1) ! extend in longitude fptr(i,j) = fptr(lbnd(1),j) enddo ! Set north pole: elseif (j==jnpole) then ! north pole if (ubnd(2) < jnpole-2) & ! need lats nlat-1:nlat for north pole call shutdown('esmf_set2d_geo need j=nlat-1:nlat') fptr(lbnd(1),j) = & (9._r8*dot_product(unitv,fpole_jpm2(:,3,1))- & dot_product(unitv,fpole_jpm2(:,4,1))) & / (8._r8*r8_nlon) do i=lbnd(1)+1,ubnd(1) ! extend in longitude fptr(i,j) = fptr(lbnd(1),j) enddo endif ! north or south pole enddo ! lat loop for poles end subroutine esmf_set2d_geo !----------------------------------------------------------------------- subroutine esmf_set3d_geo(fields,fnames,f,nf,ilev0,ilev1,& ilon0,ilon1,ilat0,ilat1) use my_mpi,only: mp_geopole_3d ! ! Set values of a 3d ESMF field on geographic source grid, prior to ! geographic to magnetic grid transformation. ! Periodic points are excluded, geographic poles are at j==jspole and jnpole ! Note dimension order changes from input (k,i,j) to output (i,j,k). ! ! Args: integer,intent(in) :: nf type(ESMF_Field) ,intent(in) :: fields(nf) ! esmf fields on geo grid character(len=*) ,intent(in) :: fnames(nf) ! field names ! ! f is input data on model subdomains (including periodic points) ! (note esmf source field excludes periodic points) ! integer,intent(in) :: ilev0,ilev1,ilon0,ilon1,ilat0,ilat1 real(r8),intent(in) :: f(ilev0:ilev1,ilon0:ilon1,ilat0:ilat1,nf) ! ! Local: integer :: i,ii,j,k,rc,n,istat integer,save :: ncalls=0 integer,save :: lbnd(3),ubnd(3) ! lower,upper bounds of 3d field ! ! fptr is esmf pointer (i,j,k) to 3d field, set by this subroutine real(ESMF_KIND_R8),pointer :: fptr(:,:,:) real(r8),allocatable,save :: ftmp(:,:,:,:) ! esmf bounds, plus nf real(r8) :: fpole_jpm2(nlon,4,ilev1,nf) ! ! This routine is called first w/ nf=1, then w/ nf=7, so could set ! mxnf=7 and allocate mxnf, but its only called twice, so am ! allocating/deallocating each time: ! call ESMF_FieldGet(fields(1),localDe=0,farrayPtr=fptr,& computationalLBound=lbnd,computationalUBound=ubnd,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(a,i4)") '>>> esmf_set3d_geo: error from ESMF_FieldGet: rc=',rc call shutdown('esmf_set3d_geo: ESMF_FieldGet') endif allocate(ftmp(lbnd(1):ubnd(1),lbnd(2):ubnd(2),lbnd(3):ubnd(3),& nf),stat=istat) if (istat /= 0) & call shutdown('esmf_set3d_ge: error allocating ftmp') ! ! Fields loop: do n=1,nf ftmp(:,:,:,n) = 0._r8 ! ! Set interior latitudes (ftmp(i,j,k,n) <- f(k,i,j,n)) ! ftmp excludes periodic points. ! do j=lbnd(2),ubnd(2) ! lat if (j /= jspole .and. j /= jnpole) then ! interior latitudes (not poles) do i=lbnd(1),ubnd(1) ! lon ii = i if (trim(model_name)=='TIMEGCM') ii = i+2 do k=lbnd(3),ubnd(3) ! lev ftmp(i,j,k,n) = f(k,ii,j,n) ! data is at i+2 (3->74) enddo ! lev enddo ! lon endif ! poles or interior enddo ! lat enddo ! n=1,nf ! ! Single call mp_geopole_3d to get global lons surrounding the 2 lats ! adjacent to each pole for each field. ! subroutine mp_geopole_3d(f,ilon0,ilon1,ilat0,ilat1,nlevs, ! nglblon,jspole,jnpole,fpole_jpm2,nf,fnames) ! Returns fpole_jpm2(nglblon,1->4,nlevs,nf) with global lons surrounding ! poles +/- 2 as follows: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) ! call mp_geopole_3d(ftmp,lbnd(1),ubnd(1),lbnd(2),ubnd(2),& nlev,nlon,jspole,jnpole,fpole_jpm2,nf,fnames) ! ! Set poles using global lons from mp_geopole_3d: do n=1,nf do j=lbnd(2),ubnd(2) ! lat ! ! Set south pole: if (j==jspole) then ! south pole if (lbnd(2) > jspole+2) & ! need lats 1,2 for south pole call shutdown('esmf_set3d_geo need j=1,2') if (trim(fnames(n))/='ADOTV1'.and.trim(fnames(n))/='ADOTV2' & .and.trim(fnames(n))/='PHIG3D') then do k=lbnd(3),ubnd(3) ! lev (1,nlev) ftmp(lbnd(1),j,k,n) = & (9._r8*dot_product(unitv,fpole_jpm2(:,1,k,n))- & dot_product(unitv,fpole_jpm2(:,2,k,n))) & / (8._r8*r8_nlon) do i=lbnd(1)+1,ubnd(1) ! extend in longitude ftmp(i,j,k,n) = ftmp(lbnd(1),j,k,n) enddo enddo ! south pole levs endif ! ADOTV or not ! ! Set north pole: elseif (j==jnpole) then ! north pole if (ubnd(2) < jnpole-2) & ! need lats nlat-1:nlat for north pole call shutdown('esmf_set3d_geo need j=nlat-1:nlat') if (trim(fnames(n))/='ADOTV1'.and.trim(fnames(n))/='ADOTV2' & .and.trim(fnames(n))/='PHIG3D') then do k=lbnd(3),ubnd(3) ! lev (1,nlev) ftmp(lbnd(1),j,k,n) = & (9._r8*dot_product(unitv,fpole_jpm2(:,3,k,n))- & dot_product(unitv,fpole_jpm2(:,4,k,n))) & / (8._r8*r8_nlon) do i=lbnd(1)+1,ubnd(1) ! extend in longitude ftmp(i,j,k,n) = ftmp(lbnd(1),j,k,n) enddo enddo ! south pole levs ! ! Set poles of poten (from source history) to match dynpot of odyn ! (sub dynpotmag in magfield.F) elseif (trim(fnames(n))=='PHIG3D') then do k=lbnd(3),ubnd(3) ! lev (1,nlev) ftmp(lbnd(1),j,k,n) = 1._r8/r8_nlon* & dot_product(unitv,fpole_jpm2(:,3,k,n)) do i=lbnd(1)+1,ubnd(1) ! extend in longitude ftmp(i,j,k,n) = ftmp(lbnd(1),j,k,n) enddo enddo else do k=lbnd(3),ubnd(3) ! lev (1,nlev) ftmp(lbnd(1),j,k,n) = & (9._r8*fpole_jpm2(1,3,k,n) - fpole_jpm2(1,4,k,n))/8._r8 do i=lbnd(1)+1,ubnd(1) ! extend in longitude ftmp(i,j,k,n) = & (9._r8*fpole_jpm2(i,3,k,n) - fpole_jpm2(i,4,k,n))/8._r8 enddo enddo ! south pole levs endif ! ADOTV or not endif ! north or south pole enddo ! j=lbnd(2),ubnd(2) ! ! Get and set pointer to the field: call ESMF_FieldGet(fields(n),localDe=0,farrayPtr=fptr, & computationalLBound=lbnd,computationalUBound=ubnd,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(a,i4)") '>>> esmf_set3d_geo: error from ESMF_FieldGet: rc=',rc call shutdown('esmf_set3d_geo:ESMF_FieldGet') endif fptr(:,:,:) = ftmp(:,:,:,n) enddo ! n=1,nf deallocate(ftmp) end subroutine esmf_set3d_geo !----------------------------------------------------------------------- subroutine esmf_set3d_mag(fields,fnames,f,nf,ilev0,ilev1,ilon0,ilon1,ilat0,ilat1) ! ! Set values of a 3d ESMF field on magnetic grid, prior to magnetic to ! geographic grid transformation. ! ! Args: integer,intent(in) :: nf type(ESMF_Field) ,intent(in) :: fields(nf) ! esmf fields on mag grid character(len=*) ,intent(in) :: fnames(nf) ! field names ! ! f is input data on model subdomains: ! integer,intent(in) :: ilev0,ilev1,ilon0,ilon1,ilat0,ilat1 real(r8),intent(in) :: f(ilon0:ilon1,ilat0:ilat1,ilev0:ilev1,nf) ! ! Local: integer :: i,j,k,rc,n,istat integer,save :: ncalls=0,lbnd(3),ubnd(3) ! lower,upper bounds of 3d field ! ! fptr is esmf pointer (i,j,k) to 3d field, set by this subroutine real(ESMF_KIND_R8),pointer :: fptr(:,:,:) ! ! Fields loop: do n=1,nf call ESMF_FieldGet(fields(n),localDe=0,farrayPtr=fptr,& computationalLBound=lbnd,computationalUBound=ubnd,rc=rc) fptr(:,:,:) = 0._r8 ! ! Set ESMF pointer: ! do j=lbnd(2),ubnd(2) ! lat do i=lbnd(1),ubnd(1) ! lon do k=lbnd(3),ubnd(3) ! lev fptr(i,j,k) = f(i,j,k,n) enddo ! mlev enddo ! mlon enddo ! mlat enddo ! n=1,nf end subroutine esmf_set3d_mag !----------------------------------------------------------------------- subroutine esmf_get_2dfield(field, fptr, name) ! ! Get pointer to 2d esmf field (i,j): ! ! Args: type(ESMF_field),intent(in) :: field real(r8),pointer,dimension(:,:),intent(out) :: fptr character(len=*),intent(in) :: name ! ! Local: integer :: rc character(len=80) :: errmsg call ESMF_FieldGet(field,localDe=0,farrayPtr=fptr,rc=rc) if (rc /= ESMF_SUCCESS) then write(errmsg,"('esmf_get_field 2d field ',a)") trim(name) call shutdown(errmsg) endif end subroutine esmf_get_2dfield !----------------------------------------------------------------------- subroutine esmf_get_3dfield(field, fptr, name) ! ! Get pointer to 3d esmf field (i,j,k): ! ! Args: type(ESMF_field),intent(in) :: field real(r8),pointer,dimension(:,:,:),intent(out) :: fptr character(len=*),intent(in) :: name ! ! Local: integer :: rc,k,lbnd(3),ubnd(3) character(len=80) :: errmsg call ESMF_FieldGet(field,localDe=0,farrayPtr=fptr, & computationalLBound=lbnd,computationalUBound=ubnd,rc=rc) if (rc /= ESMF_SUCCESS) then write(errmsg,"('esmf_get_field 3d field ',a)") trim(name) call shutdown(errmsg) endif end subroutine esmf_get_3dfield !----------------------------------------------------------------------- subroutine esmf_regrid(srcfield,dstfield,direction,ndim) ! ! Args: integer :: ndim type(ESMF_Field),intent(inout) :: srcfield,dstfield character(len=*),intent(in) :: direction ! ! Local: integer :: rc type(ESMF_RouteHandle) :: routehandle ! ! Direction is either geo2mag or mag2geo. ! Use corresponding route handle (module data) ! select case(trim(direction)) case ('geo2mag') routehandle = routehandle_geo2mag if (ndim==2) then ! ! Do sparse matrix multiply for 2d geo2mag. This replaces old ! call to ESMF_FieldRegrid). Note routehandle_geo2mag may be made ! from apex or esmf indices and weights (see call ESMF_FieldSMMStore ! in sub esmf_init) ! routehandle = routehandle_geo2mag_2d call ESMF_FieldSMM(srcfield,dstfield,routehandle,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(/,4a,i4)") '>>> esmf_regrid: error return from ',& 'ESMF_FieldSMM for 2d ',trim(direction),': rc=',rc call shutdown('ESMF_FieldRegrid') endif else ! 3d geo2mag ! ! Do sparse matrix multiply for 3d geo2mag. Note routehandle_geo2mag ! may be made from apex or esmf indices and weights (see call ! ESMF_FieldSMMStore in sub esmf_init) ! routehandle = routehandle_geo2mag call ESMF_FieldSMM(srcfield,dstfield,routehandle,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(/,4a,i4)") '>>> esmf_regrid: error return from ',& 'ESMF_FieldSMM for 3d ',trim(direction),': rc=',rc call shutdown('ESMF_FieldSMM') endif endif ! ! Do sparse matrix multiply for 3d mag2geo. Note routehandle_mag2geo ! may be made from apex or esmf indices and weights (see call ! ESMF_FieldSMMStore in sub esmf_init) ! case ('mag2geo') routehandle = routehandle_mag2geo call ESMF_FieldSMM(srcfield,dstfield,routehandle,rc=rc) if (rc /= ESMF_SUCCESS) then write(iulog,"(/,4a,i4)") '>>> esmf_regrid: error return from ',& 'ESMF_FieldSMM for 3d ',trim(direction),': rc=',rc call shutdown('ESMF_FieldSMM') endif case default write(iulog,"('>>> esmf_regrid: bad direction=',a)") trim(direction) call shutdown('esmf_regrid') end select end subroutine esmf_regrid !----------------------------------------------------------------------- end module my_esmf