Changeset 1319 for LMDZ4/trunk/libf/dyn3d
- Timestamp:
- Feb 23, 2010, 10:29:54 PM (15 years ago)
- Location:
- LMDZ4/trunk/libf/dyn3d
- Files:
-
- 3 edited
- 4 moved
-
ce0l.F90 (moved) (moved from LMDZ4/trunk/libf/dyn3d/ce0l.F) (1 diff)
-
comdissnew.h (modified) (2 diffs)
-
conf_gcm.F (modified) (2 diffs)
-
etat0_netcdf.F90 (moved) (moved from LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F) (3 diffs)
-
limit_netcdf.F90 (moved) (moved from LMDZ4/trunk/libf/dyn3d/limit_netcdf.F) (1 diff)
-
logic.h (modified) (2 diffs)
-
startvar.F90 (moved) (moved from LMDZ4/trunk/libf/dyn3d/startvar.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/trunk/libf/dyn3d/ce0l.F90
r1318 r1319 2 2 ! $Id$ 3 3 ! 4 PROGRAM create_etat0_limit 4 !------------------------------------------------------------------------------- 5 ! 6 PROGRAM ce0l 7 ! 8 !------------------------------------------------------------------------------- 9 ! Purpose: Calls etat0, creates initial states and limit_netcdf 10 ! 11 ! interbar=.T. for barycentric interpolation inter_barxy 12 ! extrap =.T. for data extrapolation, like for the SSTs when file does not 13 ! contain ocean points only. 14 ! oldice =.T. for old-style ice, obtained using grille_m (grid_atob). 15 ! masque is created in etat0, passed to limit to ensure consistancy. 16 !------------------------------------------------------------------------------- 5 17 #ifdef CPP_EARTH 6 18 ! This prog. is designed to work for Earth 7 USE dimphy 8 USE comgeomphy 9 USE infotrac 19 USE dimphy 20 USE comgeomphy 21 USE infotrac 22 10 23 #ifdef CPP_IOIPSL 11 use ioipsl, only: ioconf_calendar24 USE ioipsl, ONLY: ioconf_calendar 12 25 #endif 13 IMPLICIT NONE14 c15 c16 c Programme d'appel a etat0, creation des etats initiaux et limit_netcdf17 c18 c19 c interbar = .T . si appel a interpol. barycentrique inter_barxy20 c21 c extrap = .T . si on fait une extrapolation de donnees , comme pour22 c les SST lorsque le fichier ne contient pas uniquement des points23 c oceaniques .24 c25 c oldice = .T. si l'on veut garder les anciennes glaces , obtenues26 c par grille_m ( grid_atob ) .27 c28 c on cree le masque dans etat0 que l'on passe ensuite dans limit pour29 c garder les coherences30 26 31 LOGICAL interbar, extrap , oldice 32 PARAMETER ( interbar = .true. , extrap = .FALSE. , oldice=.false.) 27 #endif 28 IMPLICIT NONE 29 #ifndef CPP_EARTH 30 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' 31 #else 32 !------------------------------------------------------------------------------- 33 ! Local variables: 34 LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE. 33 35 #include "dimensions.h" 34 36 #include "paramet.h" 35 37 #include "indicesol.h" 36 #include "control.h" 37 REAL :: masque(iip1,jjp1) 38 ! REAL :: pctsrf(iim*(jjm-1)+2, nbsrf) 38 #include "iniprint.h" 39 #include "control.h" 40 #include "temps.h" 41 #include "logic.h" 42 INTEGER, PARAMETER :: longcles=20 43 REAL, DIMENSION(longcles) :: clesphy0 44 REAL, DIMENSION(iip1,jjp1) :: masque 45 CHARACTER(LEN=15) :: calnd 46 !------------------------------------------------------------------------------- 47 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 39 48 40 IF (config_inca /= 'none') THEN 41 #ifdef INCA 42 call init_const_lmdz( 43 $ nbtr,anneeref,dayref, 44 $ iphysiq, day_step,nday) 45 #endif 46 print *, 'nbtr =' , nbtr 47 END IF 48 49 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 50 PRINT *,'---> klon=',klon 51 call InitComgeomphy 49 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 50 WRITE(lunout,*)'---> klon=',klon 51 CALL InitComgeomphy 52 52 53 53 #ifdef CPP_IOIPSL 54 call ioconf_calendar('360d') 54 SELECT CASE(calend) 55 CASE('earth_360d');CALL ioconf_calendar('360d'); calnd='a 360 jours/an' 56 CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='a 365 jours/an' 57 CASE('earth_366d');CALL ioconf_calendar('366d'); calnd='bissextile' 58 CASE('gregorian'); CALL ioconf_calendar('gregorian'); calnd='gregorien' 59 CASE('standard'); CALL ioconf_calendar('gregorian'); calnd='gregorien' 60 CASE('julian'); CALL ioconf_calendar('julian'); calnd='julien' 61 CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian') 62 !--- DC Bof... => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian 63 CASE DEFAULT 64 CALL abort_gcm('ce0l','Mauvais choix de calendrier',1) 65 END SELECT 66 WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre '//TRIM(calnd) 55 67 #endif 56 68 57 WRITE(6,*) ' ********************* ' 58 WRITE(6,*) ' interbar = ',interbar 59 CALL etat0_netcdf ( interbar, masque ) 60 c 61 WRITE(6,1) 62 WRITE(6,*) ' ********************* ' 63 WRITE(6,*) ' *** Limit_netcdf *** ' 64 WRITE(6,*) ' ********************* ' 65 WRITE(6,1) 66 67 c 68 CALL limit_netcdf ( interbar, extrap , oldice, masque) 69 IF (config_inca /= 'none') THEN 70 #ifdef INCA 71 CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday) 72 CALL init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0) 73 WRITE(lunout,*)'nbtr =' , nbtr 74 #endif 75 END IF 69 76 70 1 FORMAT(//) 77 WRITE(lunout,'(//)') 78 WRITE(lunout,*) ' ********************* ' 79 WRITE(lunout,*) ' *** etat0_netcdf *** ' 80 WRITE(lunout,*) ' ********************* ' 81 WRITE(lunout,'(//)') 82 WRITE(lunout,*) ' interbar = ',interbar 83 CALL etat0_netcdf(interbar,masque,ok_etat0) 84 85 IF(ok_limit) THEN 86 WRITE(lunout,'(//)') 87 WRITE(lunout,*) ' ********************* ' 88 WRITE(lunout,*) ' *** Limit_netcdf *** ' 89 WRITE(lunout,*) ' ********************* ' 90 WRITE(lunout,'(//)') 91 CALL limit_netcdf(interbar,extrap,oldice,masque) 92 END IF 71 93 72 94 #endif 73 ! of #ifdef CPP_EARTH 74 STOP 75 END 95 ! of #ifndef CPP_EARTH #else 96 STOP 97 98 END PROGRAM ce0l 99 ! 100 !------------------------------------------------------------------------------- -
LMDZ4/trunk/libf/dyn3d/comdissnew.h
r524 r1319 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 c----------------------------------------------------------------------- 5 c INCLUDE comdissnew.h 4 ! 5 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre 6 ! veillez à n'utiliser que des ! pour les commentaires 7 ! et à bien positionner les & des lignes de continuation 8 ! (les placer en colonne 6 et en colonne 73) 9 ! 10 !----------------------------------------------------------------------- 11 ! INCLUDE 'comdissnew.h' 6 12 7 COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv, 8 1tetagrot,tetatemp,coefdis13 COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv, & 14 & tetagrot,tetatemp,coefdis 9 15 10 16 LOGICAL lstardis … … 12 18 REAL tetagdiv, tetagrot, tetatemp, coefdis 13 19 14 c 15 c... Les parametres de ce common comdissnew sont lues par defrun_new16 csur le fichier run.def ....17 c 18 c-----------------------------------------------------------------------20 ! 21 ! ... Les parametres de ce common comdissnew sont lues par defrun_new 22 ! sur le fichier run.def .... 23 ! 24 !----------------------------------------------------------------------- -
LMDZ4/trunk/libf/dyn3d/conf_gcm.F
r1279 r1319 769 769 CALL getin('ok_gradsfile',ok_gradsfile) 770 770 771 !Config Key = ok_limit 772 !Config Desc = creation des fichiers limit dans create_etat0_limit 773 !Config Def = y 774 !Config Help = production du fichier limit.nc requise 775 776 ok_limit = .TRUE. 777 CALL getin('ok_limit',ok_limit) 778 779 !Config Key = ok_etat0 780 !Config Desc = creation des fichiers etat0 dans create_etat0_limit 781 !Config Def = y 782 !Config Help = production des fichiers start.nc, startphy.nc requise 783 784 ok_etat0 = .TRUE. 785 CALL getin('ok_etat0',ok_etat0) 786 771 787 write(lunout,*)' #########################################' 772 write(lunout,*)' Configuration des parametres du gcm: ' 788 write(lunout,*)' Configuration des parametres de create_etat0' 789 & //'_limit: ' 773 790 write(lunout,*)' planet_type = ', planet_type 774 791 write(lunout,*)' calend = ', calend … … 809 826 write(lunout,*)' ok_strato = ', ok_strato 810 827 write(lunout,*)' ok_gradsfile = ', ok_gradsfile 828 write(lunout,*)' ok_limit = ', ok_limit 829 write(lunout,*)' ok_etat0 = ', ok_etat0 811 830 c 812 831 RETURN -
LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90
r1318 r1319 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 c 5 c 6 SUBROUTINE etat0_netcdf (interbar, masque) 7 #ifdef CPP_EARTH 8 USE startvar 9 USE ioipsl 10 USE dimphy 11 USE infotrac 12 USE fonte_neige_mod 13 USE pbl_surface_mod 14 USE phys_state_var_mod 15 USE filtreg_mod 16 use regr_lat_time_climoz_m, only: regr_lat_time_climoz 17 use conf_phys_m, only: conf_phys 4 !------------------------------------------------------------------------------- 5 ! 6 SUBROUTINE etat0_netcdf(ib, masque, letat0) 7 ! 8 !------------------------------------------------------------------------------- 9 ! Purpose: Creates initial states 10 !------------------------------------------------------------------------------- 11 ! Note: This routine is designed to work for Earth 12 !------------------------------------------------------------------------------- 13 #ifdef CPP_EARTH 14 USE startvar 15 USE ioipsl 16 USE dimphy 17 USE infotrac 18 USE fonte_neige_mod 19 USE pbl_surface_mod 20 USE phys_state_var_mod 21 USE filtreg_mod 22 USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz 23 USE conf_phys_m, ONLY: conf_phys 24 USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR 18 25 #endif 19 !#endif of #ifdef CPP_EARTH 20 use netcdf, only: nf90_open, NF90_NOWRITE, nf90_close 21 ! 22 IMPLICIT NONE 23 ! 26 IMPLICIT NONE 27 !------------------------------------------------------------------------------- 28 ! Arguments: 24 29 #include "dimensions.h" 25 30 #include "paramet.h" 26 ! 27 ! 28 ! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, 29 ! .KLON=KFDIA-KIDIA+1,KLEV=llm 30 ! 31 #ifdef CPP_EARTH 31 #include "iniprint.h" 32 LOGICAL, INTENT(IN) :: ib ! barycentric interpolat. 33 REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask 34 LOGICAL, INTENT(IN) :: letat0 ! F: masque only required 35 #ifndef CPP_EARTH 36 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' 37 #else 38 !------------------------------------------------------------------------------- 39 ! Local variables: 32 40 #include "comgeom2.h" 33 41 #include "comvert.h" … … 36 44 #include "dimsoil.h" 37 45 #include "temps.h" 38 #endif 39 !#endif of #ifdef CPP_EARTH 40 ! arguments: 41 LOGICAL interbar 42 REAL :: masque(iip1,jjp1) 43 44 #ifdef CPP_EARTH 45 ! local variables: 46 REAL :: latfi(klon), lonfi(klon) 47 REAL :: orog(iip1,jjp1), rugo(iip1,jjp1) 48 REAL :: psol(iip1, jjp1), phis(iip1, jjp1) 49 REAL :: p3d(iip1, jjp1, llm+1) 50 REAL :: uvent(iip1, jjp1, llm) 51 REAL :: vvent(iip1, jjm, llm) 52 REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm) 53 REAL :: qsat(iip1, jjp1, llm) 54 REAL,ALLOCATABLE :: q3d(:, :, :,:) 55 REAL :: tsol(klon), qsol(klon), sn(klon) 56 !! REAL :: tsolsrf(klon,nbsrf) 57 real qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 58 REAL :: albe(klon,nbsrf), evap(klon,nbsrf) 59 REAL :: alblw(klon,nbsrf) 60 REAL :: tsoil(klon,nsoilmx,nbsrf) 61 REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf) 62 REAL :: rugmer(klon) 63 REAL :: qd(iip1, jjp1, llm) 64 REAL :: run_off_lic_0(klon) 65 ! declarations pour lecture glace de mer 66 REAL :: rugv(klon) 67 INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret 68 INTEGER :: itaul(1), fid 69 REAL :: lev(1), date 70 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic 71 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic 72 REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic 73 REAL :: flic_tmp(iip1, jjp1) 74 REAL :: champint(iim, jjp1) 75 ! 76 77 CHARACTER(len=80) :: varname 78 ! 79 INTEGER :: i,j, ig, l, ji,ii1,ii2 80 REAL :: xpi 81 ! 82 REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) 83 REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1) 84 REAL :: workvar(iip1,jjp1,llm) 85 ! 86 REAL :: prefkap, unskap 87 ! 88 real :: time_step,t_ops,t_wrt 46 REAL, DIMENSION(klon) :: tsol, qsol 47 REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0 48 REAL, DIMENSION(iip1,jjp1) :: orog, rugo, psol, phis 49 REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d 50 REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd 51 REAL, DIMENSION(iip1,jjm ,llm) :: vvent 52 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: q3d 53 REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf, evap 54 REAL, DIMENSION(klon,nbsrf) :: frugs, agesno 55 REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil 56 57 !--- Local variables for sea-ice reading: 58 INTEGER :: iml_lic, jml_lic, llm_tmp 59 INTEGER :: ttm_tmp, iret, fid 60 INTEGER, DIMENSION(1) :: itaul 61 REAL, DIMENSION(1) :: lev 62 REAL :: date 63 REAL, DIMENSION(:,:), ALLOCATABLE :: lon_lic, lat_lic, fraclic 64 REAL, DIMENSION(:), ALLOCATABLE :: dlon_lic, dlat_lic 65 REAL, DIMENSION(iip1,jjp1) :: flic_tmp 66 67 !--- Misc 68 CHARACTER(LEN=80) :: x, fmt 69 INTEGER :: i, j, l, ji 70 REAL, DIMENSION(iip1,jjp1,llm) :: alpha, beta, pk, pls, y 71 REAL, DIMENSION(ip1jmp1) :: pks 89 72 90 73 #include "comdissnew.h" … … 93 76 #include "clesphys.h" 94 77 95 INTEGER :: longcles 96 PARAMETER ( longcles = 20 ) 97 REAL :: clesphy0 ( longcles ) 98 REAL :: p(iip1,jjp1,llm) 99 INTEGER :: itau, iday 100 REAL :: masse(iip1,jjp1,llm) 101 REAL :: xpn,xps,xppn(iim),xpps(iim) 102 real :: time 103 REAL :: phi(ip1jmp1,llm) 104 REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 105 REAL :: w(ip1jmp1,llm) 106 REAL ::phystep 107 CC REAL :: rugsrel(iip1*jjp1) 108 REAL :: fder(klon) 109 !! real zrel(iip1*jjp1),chmin,chmax 110 111 !! CHARACTER(len=80) :: visu_file 112 INTEGER :: visuid 113 114 ! pour la lecture du fichier masque ocean 115 integer :: nid_o2a 116 logical :: couple = .false. 117 INTEGER :: iml_omask, jml_omask 118 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask 119 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask 120 REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp 121 real, dimension(klon) :: ocemask_fi 122 integer :: isst(klon-2) 123 real zx_tmp_2d(iim,jjp1) 124 125 REAL :: dummy 126 127 logical :: ok_newmicro 128 integer :: iflag_radia 129 logical :: ok_journe, ok_mensuel, ok_instan, ok_hf 130 logical :: ok_LES 131 LOGICAL :: ok_ade, ok_aie, aerosol_couple, new_aod 132 INTEGER :: flag_aerosol 133 REAL :: bl95_b0, bl95_b1 134 real :: fact_cldcon, facttemps,ratqsbas,ratqshaut 135 real :: tau_ratqs 136 integer :: iflag_cldcon 137 integer :: iflag_ratqs 138 integer :: iflag_coupl 139 integer :: iflag_clos 140 integer :: iflag_wake 141 integer :: iflag_thermals,nsplit_thermals 142 real :: tau_thermals 143 integer :: iflag_thermals_ed,iflag_thermals_optflux 144 REAL :: solarlong0 145 real :: seuil_inversion 146 147 integer read_climoz ! read ozone climatology 148 C Allowed values are 0, 1 and 2 149 C 0: do not read an ozone climatology 150 C 1: read a single ozone climatology that will be used day and night 151 C 2: read two ozone climatologies, the average day and night 152 C climatology and the daylight climatology 153 154 ! 155 ! Constantes 156 ! 157 pi = 4. * ATAN(1.) 158 rad = 6371229. 159 omeg = 4.* ASIN(1.)/(24.*3600.) 160 g = 9.8 161 daysec = 86400. 162 kappa = 0.2857143 163 cpp = 1004.70885 164 ! 165 preff = 101325. 166 pa = 50000. 167 unskap = 1./kappa 168 ! 169 jmp1 = jjm + 1 170 ! 171 ! Construct a grid 172 ! 173 174 ! CALL defrun_new(99,.TRUE.,clesphy0) 175 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 176 call conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & 177 & solarlong0,seuil_inversion, & 178 & fact_cldcon, facttemps,ok_newmicro,iflag_radia, & 179 & iflag_cldcon, & 180 & iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & 181 & ok_ade, ok_aie, aerosol_couple, & 182 & flag_aerosol, new_aod, & 183 & bl95_b0, bl95_b1, & 184 & iflag_thermals,nsplit_thermals,tau_thermals, & 185 & iflag_thermals_ed,iflag_thermals_optflux, & 186 & iflag_coupl,iflag_clos,iflag_wake, read_climoz ) 78 REAL, DIMENSION(iip1,jjp1,llm) :: masse 79 INTEGER :: itau, iday 80 REAL :: xpn, xps, time, phystep 81 REAL, DIMENSION(iim) :: xppn, xpps 82 REAL, DIMENSION(ip1jmp1,llm) :: pbaru, phi, w 83 REAL, DIMENSION(ip1jm ,llm) :: pbarv 84 REAL, DIMENSION(klon) :: fder 85 86 !--- Local variables for ocean mask reading: 87 INTEGER :: nid_o2a, iml_omask, jml_omask 88 LOGICAL :: couple=.FALSE. 89 REAL, DIMENSION(:,:), ALLOCATABLE :: lon_omask, lat_omask, ocemask, ocetmp 90 REAL, DIMENSION(:), ALLOCATABLE :: dlon_omask,dlat_omask 91 REAL, DIMENSION(klon) :: ocemask_fi 92 INTEGER, DIMENSION(klon-2) :: isst 93 REAL, DIMENSION(iim,jjp1) :: zx_tmp_2d 94 REAL :: dummy 95 LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf 96 LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod 97 INTEGER :: iflag_radia, flag_aerosol 98 REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut 99 REAL :: tau_ratqs 100 INTEGER :: iflag_cldcon, iflag_ratqs, iflag_coupl, iflag_clos, iflag_wake 101 INTEGER :: iflag_thermals, nsplit_thermals 102 INTEGER :: iflag_thermals_ed, iflag_thermals_optflux 103 REAL :: tau_thermals, solarlong0, seuil_inversion 104 INTEGER :: read_climoz ! read ozone climatology 105 ! Allowed values are 0, 1 and 2 106 ! 0: do not read an ozone climatology 107 ! 1: read a single ozone climatology that will be used day and night 108 ! 2: read two ozone climatologies, the average day and night 109 ! climatology and the daylight climatology 110 !------------------------------------------------------------------------------- 111 !--- Constants 112 pi = 4. * ATAN(1.) 113 rad = 6371229. 114 daysec = 86400. 115 omeg = 2.*pi/daysec 116 g = 9.8 117 kappa = 0.2857143 118 cpp = 1004.70885 119 preff = 101325. 120 pa = 50000. 121 jmp1 = jjm + 1 122 123 !--- CONSTRUCT A GRID 124 CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & 125 solarlong0,seuil_inversion, & 126 fact_cldcon, facttemps,ok_newmicro,iflag_radia, & 127 iflag_cldcon, & 128 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & 129 ok_ade, ok_aie, aerosol_couple, & 130 flag_aerosol, new_aod, & 131 bl95_b0, bl95_b1, & 132 iflag_thermals,nsplit_thermals,tau_thermals, & 133 iflag_thermals_ed,iflag_thermals_optflux, & 134 iflag_coupl,iflag_clos,iflag_wake, read_climoz ) 187 135 188 136 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value) 189 co2_ppm0 = co2_ppm 190 191 dtvr = daysec/FLOAT(day_step) 192 print*,'dtvr',dtvr 193 194 CALL iniconst() 195 CALL inigeom() 196 197 ! Initialisation pour traceurs 198 call infotrac_init 199 ALLOCATE(q3d(iip1, jjp1, llm, nqtot)) 200 201 CALL inifilr() 202 CALL phys_state_var_init(read_climoz) 203 ! 204 latfi(1) = ASIN(1.0) 205 DO j = 2, jjm 206 DO i = 1, iim 207 latfi((j-2)*iim+1+i)= rlatu(j) 208 ENDDO 209 ENDDO 210 latfi(klon) = - ASIN(1.0) 211 ! 212 lonfi(1) = 0.0 213 DO j = 2, jjm 214 DO i = 1, iim 215 lonfi((j-2)*iim+1+i) = rlonv(i) 216 ENDDO 217 ENDDO 218 lonfi(klon) = 0.0 219 ! 220 xpi = 2.0 * ASIN(1.0) 221 DO ig = 1, klon 222 latfi(ig) = latfi(ig) * 180.0 / xpi 223 lonfi(ig) = lonfi(ig) * 180.0 / xpi 224 ENDDO 225 ! 226 rlat(1) = ASIN(1.0) 227 DO j = 2, jjm 228 DO i = 1, iim 229 rlat((j-2)*iim+1+i)= rlatu(j) 230 ENDDO 231 ENDDO 232 rlat(klon) = - ASIN(1.0) 233 ! 234 rlon(1) = 0.0 235 DO j = 2, jjm 236 DO i = 1, iim 237 rlon((j-2)*iim+1+i) = rlonv(i) 238 ENDDO 239 ENDDO 240 rlon(klon) = 0.0 241 ! 242 xpi = 2.0 * ASIN(1.0) 243 DO ig = 1, klon 244 rlat(ig) = rlat(ig) * 180.0 / xpi 245 rlon(ig) = rlon(ig) * 180.0 / xpi 246 ENDDO 247 ! 248 249 250 251 C 252 C En cas de simulation couplee, lecture du masque ocean issu du modele ocean 253 C utilise pour calculer les poids et pour assurer l'adequation entre les 254 C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque 255 C a partir du fichier relief 256 C 257 258 write(*,*)'Essai de lecture masque ocean' 259 iret = nf90_open("o2a.nc", NF90_NOWRITE, nid_o2a) 260 if (iret .ne. 0) then 261 write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve' 262 write(*,*)'Run force' 263 varname = 'masque' 264 masque(:,:) = 0.0 265 CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0, 266 , jjm ,rlonu,rlatv , interbar ) 267 WRITE(*,*) 'MASQUE construit : Masque' 268 WRITE(*,'(97I1)') nINT(masque(:,:)) 269 call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq) 270 WHERE (zmasq(1 : klon) .LT. EPSFRA) 271 zmasq(1 : klon) = 0. 272 END WHERE 273 WHERE (1. - zmasq(1 : klon) .LT. EPSFRA) 274 zmasq(1 : klon) = 1. 275 END WHERE 276 else 277 couple = .true. 278 iret = nf90_close(nid_o2a) 279 call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp 280 $ , nid_o2a) 281 if (iml_omask /= iim .or. jml_omask /= jjp1) then 282 write(*,*)'Dimensions non compatibles pour masque ocean' 283 write(*,*)'iim = ',iim,' iml_omask = ',iml_omask 284 write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask 285 stop 286 endif 287 ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret) 288 ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret) 289 ALLOCATE(dlon_omask(iml_omask), stat=iret) 290 ALLOCATE(dlat_omask(jml_omask), stat=iret) 291 ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret) 292 ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret) 293 CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp 294 $ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 295 CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, 296 $ ttm_tmp, 1, 1, ocetmp) 297 CALL flinclo(fid) 298 dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1) 299 dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask) 300 ocemask = ocetmp 301 if (dlat_omask(1) < dlat_omask(jml_omask)) then 302 do j = 1, jml_omask 303 ocemask(:,j) = ocetmp(:,jml_omask-j+1) 304 enddo 305 endif 306 C 307 C passage masque ocean a la grille physique 308 C 309 write(*,*)'ocemask ' 310 write(*,'(96i1)')int(ocemask) 311 ocemask_fi(1) = ocemask(1,1) 312 do j = 2, jjm 313 do i = 1, iim 314 ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j) 315 enddo 316 enddo 317 ocemask_fi(klon) = ocemask(1,jjp1) 318 zmasq = 1. - ocemask_fi 319 endif 320 321 call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 322 323 varname = 'relief' 324 ! This line needs to be replaced by a call to restget to get the values in the restart file 325 orog(:,:) = 0.0 326 CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 , 327 , jjm ,rlonu,rlatv , interbar, masque ) 328 ! 329 WRITE(*,*) 'OUT OF GET VARIABLE : Relief' 330 ! WRITE(*,'(49I1)') INT(orog(:,:)) 331 ! 332 varname = 'rugosite' 333 ! This line needs to be replaced by a call to restget to get the values in the restart file 334 rugo(:,:) = 0.0 335 CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 , 336 , jjm, rlonu,rlatv , interbar ) 337 ! 338 WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 339 ! WRITE(*,'(49I1)') INT(rugo(:,:)*10) 340 ! 341 C 342 C on initialise les sous surfaces 343 C 344 pctsrf=0. 345 c 346 varname = 'psol' 347 psol(:,:) = 0.0 348 CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 , 349 , jjm ,rlonu,rlatv , interbar ) 350 ! 351 ! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM 352 ! anyway. 353 ! 354 ! WRITE(*,*) 'PSOL :', psol(10,20) 355 ! WRITE(*,*) ap(:), bp(:) 356 CALL pression(ip1jmp1, ap, bp, psol, p3d) 357 ! WRITE(*,*) 'P3D :', p3d(10,20,:) 358 CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar) 359 ! WRITE(*,*) 'PK:', pk(10,20,:) 360 ! 361 ! 362 ! 363 prefkap = preff ** kappa 364 ! WRITE(*,*) 'unskap, cpp, preff :', unskap, cpp, preff 365 DO l = 1, llm 366 DO j=1,jjp1 367 DO i =1, iip1 368 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap 369 ENDDO 370 ENDDO 371 ENDDO 372 ! 373 ! WRITE(*,*) 'PLS :', pls(10,20,:) 374 ! 375 varname = 'surfgeo' 376 phis(:,:) = 0.0 377 CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 , 378 , jjm ,rlonu,rlatv, interbar ) 379 ! 380 varname = 'u' 381 uvent(:,:,:) = 0.0 382 CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls, 383 . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar ) 384 ! 385 varname = 'v' 386 vvent(:,:,:) = 0.0 387 CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls, 388 . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar ) 389 ! 390 varname = 't' 391 t3d(:,:,:) = 0.0 392 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, 393 . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar ) 394 ! 395 WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), 396 . maxval(t3d(:,:,:)) 397 varname = 'tpot' 398 tpot(:,:,:) = 0.0 399 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, 400 . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar ) 401 ! 402 WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), 403 . maxval(t3d(:,:,:)) 404 WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)), 405 . maxval(pls(:,:,:)) 406 407 c Calcul de l'humidite a saturation 408 print*,'avant q_sat' 409 call q_sat(llm*jjp1*iip1,t3d,pls,qsat) 410 print*,'apres q_sat' 411 412 WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), 413 . maxval(qsat(:,:,:)) 414 ! 415 CC WRITE(*,*) 'QSAT :', qsat(10,20,:) 416 ! 417 varname = 'q' 418 qd(:,:,:) = 0.0 419 q3d(:,:,:,:) = 0.0 420 WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), 421 . maxval(qsat(:,:,:)) 422 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, 423 . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar ) 424 q3d(:,:,:,1) = qd(:,:,:) 425 ! 426 427 ! Ozone climatology: 428 if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz) 429 430 varname = 'tsol' 431 ! This line needs to be replaced by a call to restget to get the values in the restart file 432 tsol(:) = 0.0 433 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0, 434 . jjm, rlonu, rlatv , interbar ) 435 ! 436 WRITE(*,*) 'TSOL construit :' 437 ! WRITE(*,'(48I3)') INT(TSOL(2:klon)-273) 438 ! 439 varname = 'qsol' 440 qsol(:) = 0.0 441 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0, 442 . jjm, rlonu, rlatv , interbar ) 443 ! 444 varname = 'snow' 445 sn(:) = 0.0 446 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0, 447 . jjm, rlonu, rlatv , interbar ) 448 ! 449 varname = 'rads' 450 radsol(:) = 0.0 451 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0, 452 . jjm, rlonu, rlatv , interbar ) 453 ! 454 varname = 'rugmer' 455 rugmer(:) = 0.0 456 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0, 457 . jjm, rlonu, rlatv , interbar ) 458 ! 459 ! varname = 'agesno' 460 ! agesno(:) = 0.0 461 ! CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0, 462 ! . jjm, rlonu, rlatv , interbar ) 463 464 varname = 'zmea' 465 zmea(:) = 0.0 466 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0, 467 . jjm, rlonu, rlatv , interbar ) 468 469 varname = 'zstd' 470 zstd(:) = 0.0 471 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0, 472 . jjm, rlonu, rlatv , interbar ) 473 varname = 'zsig' 474 zsig(:) = 0.0 475 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0, 476 . jjm, rlonu, rlatv , interbar ) 477 varname = 'zgam' 478 zgam(:) = 0.0 479 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0, 480 . jjm, rlonu, rlatv , interbar ) 481 varname = 'zthe' 482 zthe(:) = 0.0 483 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0, 484 . jjm, rlonu, rlatv , interbar ) 485 varname = 'zpic' 486 zpic(:) = 0.0 487 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0, 488 . jjm, rlonu, rlatv , interbar ) 489 varname = 'zval' 490 zval(:) = 0.0 491 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0, 492 . jjm, rlonu, rlatv , interbar ) 493 c 494 cc rugsrel(:) = 0.0 495 cc IF(ok_orodr) THEN 496 cc DO i = 1, iip1* jjp1 497 cc rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. ) 498 cc ENDDO 499 cc ENDIF 500 501 502 C 503 C lecture du fichier glace de terre pour fixer la fraction de terre 504 C et de glace de terre 505 C 506 CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp 507 $ , fid) 508 ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret) 509 ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret) 510 ALLOCATE(dlon_lic(iml_lic), stat=iret) 511 ALLOCATE(dlat_lic(jml_lic), stat=iret) 512 ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret) 513 CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp 514 $ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) 515 CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp 516 $ , 1, 1, fraclic) 517 CALL flinclo(fid) 518 C 519 C interpolation sur la grille T du modele 520 C 521 WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 522 $ iml_lic, jml_lic 523 c 524 C sil les coordonnees sont en degres, on les transforme 525 C 526 IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN 527 lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180. 528 ENDIF 529 IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 530 lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180. 531 ENDIF 532 533 dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1) 534 dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 535 C 536 CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic 537 $ ,iim, jjp1, 538 $ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1)) 539 cx$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1) 540 flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1) 541 C 542 C passage sur la grille physique 543 C 544 CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, 545 $ pctsrf(1:klon, is_lic)) 546 C adequation avec le maque terre/mer 547 c zmasq(157) = 0. 548 WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 549 pctsrf(1 : klon, is_lic) = 0. 550 END WHERE 551 WHERE (zmasq( 1 : klon) .LT. EPSFRA) 552 pctsrf(1 : klon, is_lic) = 0. 553 END WHERE 554 pctsrf(1 : klon, is_ter) = zmasq(1 : klon) 555 DO ji = 1, klon 556 IF (zmasq(ji) .GT. EPSFRA) THEN 557 IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN 558 pctsrf(ji, is_lic) = zmasq(ji) 559 pctsrf(ji, is_ter) = 0. 560 ELSE 561 pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic) 562 IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN 563 pctsrf(ji,is_ter) = 0. 564 pctsrf(ji, is_lic) = zmasq(ji) 565 ENDIF 566 ENDIF 567 ENDIF 568 END DO 569 C 570 C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0) 571 C 572 pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon)) 573 574 575 WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA) 576 pctsrf(1 : klon, is_oce) = 0. 577 END WHERE 578 579 if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon) 580 581 isst = 0 582 where (pctsrf(2:klon-1,is_oce) >0.) isst = 1 583 C 584 C verif que somme des sous surface = 1 585 C 586 ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0) 587 $ .GT. EPSFRA) 588 IF (ji .NE. 0) THEN 589 WRITE(*,*) 'pb repartition sous maille pour ',ji,' points' 590 ENDIF 591 592 ! where (pctsrf(1:klon, is_ter) >= .5) 593 ! pctsrf(1:klon, is_ter) = 1. 594 ! pctsrf(1:klon, is_oce) = 0. 595 ! pctsrf(1:klon, is_sic) = 0. 596 ! pctsrf(1:klon, is_lic) = 0. 597 ! zmasq = 1. 598 ! endwhere 599 ! where (pctsrf(1:klon, is_lic) >= .5) 600 ! pctsrf(1:klon, is_ter) = 0. 601 ! pctsrf(1:klon, is_oce) = 0. 602 ! pctsrf(1:klon, is_sic) = 0. 603 ! pctsrf(1:klon, is_lic) = 1. 604 ! zmasq = 1. 605 ! endwhere 606 ! where (pctsrf(1:klon, is_oce) >= .5) 607 ! pctsrf(1:klon, is_ter) = 0. 608 ! pctsrf(1:klon, is_oce) = 1. 609 ! pctsrf(1:klon, is_sic) = 0. 610 ! pctsrf(1:klon, is_lic) = 0. 611 ! zmasq = 0. 612 ! endwhere 613 ! where (pctsrf(1:klon, is_sic) >= .5) 614 ! pctsrf(1:klon, is_ter) = 0. 615 ! pctsrf(1:klon, is_oce) = 0. 616 ! pctsrf(1:klon, is_sic) = 1. 617 ! pctsrf(1:klon, is_lic) = 0. 618 ! zmasq = 0. 619 ! endwhere 620 ! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 621 C 622 C verif que somme des sous surface = 1 623 C 624 ! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 ) 625 ! $ .GT. EPSFRA) 626 ! IF (ji .NE. 0) THEN 627 ! WRITE(*,*) 'pb repartition sous maille pour ',ji,' points' 628 ! ENDIF 629 630 CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d) 631 write(*,*)'zmasq = ' 632 write(*,'(96i1)')nint(zx_tmp_2d) 633 call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 634 WRITE(*,*) 'MASQUE construit : Masque' 635 WRITE(*,'(97I1)') nINT(masque(:,:)) 636 637 638 639 C Calcul intermediaire 640 c 641 CALL massdair( p3d, masse ) 642 c 643 644 print *,' ALPHAX ',alphax 645 646 DO l = 1, llm 647 DO i = 1, iim 648 xppn(i) = aire( i, 1 ) * masse( i , 1 , l ) 649 xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l ) 650 ENDDO 651 xpn = SUM(xppn)/apoln 652 xps = SUM(xpps)/apols 653 DO i = 1, iip1 654 masse( i , 1 , l ) = xpn 655 masse( i , jjp1 , l ) = xps 656 ENDDO 657 ENDDO 658 q3d(iip1,:,:,:) = q3d(1,:,:,:) 659 phis(iip1,:) = phis(1,:) 660 661 C Ecriture 662 CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , 663 * tetagdiv, tetagrot , tetatemp ) 664 print*,'sortie inidissip' 665 itau = 0 666 itau_dyn = 0 667 itau_phy = 0 668 iday = dayref +itau/day_step 669 time = real(itau-(iday-dayref)*day_step)/day_step 670 c 671 IF(time.GT.1) THEN 672 time = time - 1 673 iday = iday + 1 674 ENDIF 675 day_ref = dayref 676 annee_ref = anneeref 677 678 CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi ) 679 print*,'sortie geopot' 680 681 CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis , 682 * phi,w, pbaru,pbarv,time+iday-dayref ) 683 print*,'sortie caldyn0' 684 CALL dynredem0("start.nc",dayref,phis) 685 print*,'sortie dynredem0' 686 CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,masse , 687 . psol) 688 print*,'sortie dynredem1' 689 C 690 C Ecriture etat initial physique 691 C 692 write(*,*)'phystep ',dtvr,iphysiq,nbapp_rad 693 phystep = dtvr * FLOAT(iphysiq) 694 radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) 695 write(*,*)'phystep =', phystep, radpas 696 cIM : lecture de co2_ppm & solaire ds physiq.def 697 c co2_ppm = 348.0 698 c solaire = 1365.0 699 700 c 701 c Initialisation 702 c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs 703 c 704 ftsol(:,is_ter) = tsol 705 ftsol(:,is_lic) = tsol 706 ftsol(:,is_oce) = tsol 707 ftsol(:,is_sic) = tsol 708 snsrf(:,is_ter) = sn 709 snsrf(:,is_lic) = sn 710 snsrf(:,is_oce) = sn 711 snsrf(:,is_sic) = sn 712 falb1(:,is_ter) = 0.08 713 falb1(:,is_lic) = 0.6 714 falb1(:,is_oce) = 0.5 715 falb1(:,is_sic) = 0.6 716 falb2 = falb1 717 evap(:,:) = 0. 718 qsolsrf(:,is_ter) = 150 719 qsolsrf(:,is_lic) = 150 720 qsolsrf(:,is_oce) = 150. 721 qsolsrf(:,is_sic) = 150. 722 do i = 1, nbsrf 723 do j = 1, nsoilmx 724 tsoil(:,j,i) = tsol 725 enddo 726 enddo 727 rain_fall = 0.; snow_fall = 0. 728 solsw = 165. 729 sollw = -53. 730 t_ancien = 273.15 731 q_ancien = 0. 732 agesno = 0. 733 c 734 frugs(1:klon,is_oce) = rugmer(1:klon) 735 frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0) 736 frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0) 737 frugs(1:klon,is_sic) = 0.001 738 fder = 0.0 739 clwcon = 0.0 740 rnebcon = 0.0 741 ratqs = 0.0 742 run_off_lic_0 = 0.0 743 rugoro = 0.0 744 745 c 746 c Avant l'appel a phyredem, on initialize les modules de surface 747 c avec les valeurs qui vont etre ecrit dans startphy.nc 748 c 749 dummy = 1.0 750 pbl_tke(:,:,:) = 1.e-8 751 zmax0(:) = 40. 752 f0(:) = 1.e-5 753 ema_work1(:,:) = 0. 754 ema_work2(:,:) = 0. 755 wake_deltat(:,:) = 0. 756 wake_deltaq(:,:) = 0. 757 wake_s(:) = 0. 758 wake_cstar(:) = 0. 759 wake_fip(:) = 0. 760 761 call fonte_neige_init(run_off_lic_0) 762 call pbl_surface_init(qsol, fder, snsrf, qsolsrf, 763 $ evap, frugs, agesno, tsoil) 764 765 call phyredem("startphy.nc") 766 767 768 769 C Sortie Visu pour les champs dynamiques 770 cc if (1.eq.0 ) then 771 cc print*,'sortie visu' 772 cc time_step = 1. 773 cc t_ops = 2. 774 cc t_wrt = 2. 775 cc itau = 2. 776 cc visu_file='Etat0_visu.nc' 777 cc CALL initdynav(visu_file,dayref,anneeref,time_step, 778 cc . t_ops, t_wrt, visuid) 779 cc CALL writedynav(visuid, itau,vvent , 780 cc . uvent,tpot,pk,phi,q3d,masse,psol,phis) 781 cc else 782 print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' 783 cc endif 784 print*,'entree histclo' 785 CALL histclo 137 co2_ppm0 = co2_ppm 138 139 dtvr = daysec/FLOAT(day_step) 140 WRITE(lunout,*)'dtvr',dtvr 141 142 CALL iniconst() 143 CALL inigeom() 144 145 !--- Initializations for tracers 146 CALL infotrac_init 147 ALLOCATE(q3d(iip1,jjp1,llm,nqtot)) 148 149 CALL inifilr() 150 CALL phys_state_var_init(read_climoz) 151 152 rlat(1) = ASIN(1.0) 153 DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j); END DO 154 rlat(klon) = - ASIN(1.0) 155 rlat(:)=rlat(:)*(180.0/pi) 156 157 rlon(1) = 0.0 158 DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO 159 rlon(klon) = 0.0 160 rlon(:)=rlon(:)*(180.0/pi) 161 162 ! For a coupled simulation, the ocean mask from ocean model is used to compute 163 ! the weights an to insure ocean fractions are the same for atmosphere and ocean 164 ! Otherwise, mask is created using Relief file. 165 166 WRITE(lunout,*)'Essai de lecture masque ocean' 167 iret = NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a) 168 IF(iret/=NF90_NOERR) THEN 169 WRITE(lunout,*)'ATTENTION!! pas de fichier o2a.nc trouve' 170 WRITE(lunout,*)'Run force' 171 x='masque' 172 masque(:,:)=0.0 173 CALL startget(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, rlonu, & 174 rlatv, ib) 175 WRITE(lunout,*)'MASQUE construit : Masque' 176 WRITE(lunout,'(97I1)') nINT(masque) 177 CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq) 178 WHERE( zmasq(:)<EPSFRA) zmasq(:)=0. 179 WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1. 180 ELSE 181 couple=.true. 182 iret=NF90_CLOSE(nid_o2a) 183 CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a) 184 IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN 185 WRITE(lunout,*)'Dimensions non compatibles pour masque ocean' 186 WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask 187 WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask 188 CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc& 189 &ean',1) 190 END IF 191 ALLOCATE( ocemask(iml_omask,jml_omask), ocetmp(iml_omask,jml_omask)) 192 ALLOCATE(lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask)) 193 ALLOCATE(dlon_omask(iml_omask),dlat_omask(jml_omask)) 194 CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, & 195 lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 196 CALL flinget(fid,'OceMask',iml_omask,jml_omask,llm_tmp,ttm_tmp,1,1,ocetmp) 197 CALL flinclo(fid) 198 dlon_omask(:)=lon_omask(:,1) 199 dlat_omask(:)=lat_omask(1,:) 200 ocemask=ocetmp 201 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 202 DO j=1,jml_omask 203 ocemask(:,j) = ocetmp(:,jml_omask-j+1) 204 END DO 205 END IF 206 ALLOCATE( ocemask(iml_omask,jml_omask), ocetmp(iml_omask,jml_omask)) 207 ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask)) 208 ALLOCATE(dlon_omask(iml_omask), dlat_omask(jml_omask)) 209 CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask, & 210 lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 211 CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp, & 212 1, 1, ocetmp) 213 CALL flinclo(fid) 214 dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1) 215 dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask) 216 ocemask = ocetmp 217 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 218 DO j=1,jml_omask 219 ocemask(:,j) = ocetmp(:,jml_omask-j+1) 220 END DO 221 END IF 222 ! 223 ! Ocean mask to physical grid 224 !******************************************************************************* 225 WRITE(lunout,*)'ocemask ' 226 WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt) 227 WRITE(lunout,fmt)int(ocemask) 228 ocemask_fi(1)=ocemask(1,1) 229 DO j=2,jjm; ocemask_fi((j-2)*iim+1:iim+1)=ocemask(1:iim,j); END DO 230 ocemask_fi(klon)=ocemask(1,jjp1) 231 zmasq=1.-ocemask_fi 232 END IF 233 234 CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque) 235 236 ! The startget calls need to be replaced by a call to restget to get the 237 ! values in the restart file 238 x = 'relief'; orog(:,:) = 0.0 239 CALL startget(x,iip1,jjp1,rlonv,rlatu, orog, 0.0,jjm,rlonu,rlatv,ib,& 240 masque) 241 242 x = 'rugosite'; rugo(:,:) = 0.0 243 CALL startget(x,iip1,jjp1,rlonv,rlatu, rugo, 0.0,jjm, rlonu,rlatv,ib) 244 ! WRITE(lunout,'(49I1)') INT(orog(:,:)*10) 245 ! WRITE(lunout,'(49I1)') INT(rugo(:,:)*10) 246 247 ! Sub-surfaces initialization 248 !******************************************************************************* 249 pctsrf=0. 250 x = 'psol'; psol(:,:) = 0.0 251 CALL startget(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib) 252 ! WRITE(lunout,*) 'PSOL :', psol(10,20) 253 ! WRITE(lunout,*) ap(:), bp(:) 254 255 ! Mid-levels pressure computation 256 !******************************************************************************* 257 CALL pression(ip1jmp1, ap, bp, psol, p3d) 258 CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y) 259 pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa) 260 ! WRITE(lunout,*) 'P3D :', p3d(10,20,:) 261 ! WRITE(lunout,*) 'PK:', pk(10,20,:) 262 ! WRITE(lunout,*) 'PLS :', pls(10,20,:) 263 264 x = 'surfgeo'; phis(:,:) = 0.0 265 CALL startget(x,iip1,jjp1,rlonv,rlatu, phis, 0.0,jjm, rlonu,rlatv,ib) 266 267 x = 'u'; uvent(:,:,:) = 0.0 268 CALL startget(x,iip1,jjp1,rlonu,rlatu,llm,pls,y,uvent,0.0,jjm, rlonv,rlatv,ib) 269 270 x = 'v'; vvent(:,:,:) = 0.0 271 CALL startget(x,iip1,jjm, rlonv,rlatv,llm,pls,y,vvent,0.0,jjp1,rlonu,rlatu,ib) 272 273 x = 't'; t3d(:,:,:) = 0.0 274 CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,y,t3d, 0.0,jjm, rlonu,rlatv,ib) 275 276 x = 'tpot'; tpot(:,:,:) = 0.0 277 CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,pk,tpot,0.0,jjm, rlonu,rlatv,ib) 278 279 WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:)) 280 WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:)) 281 282 ! Humidity at saturation computation 283 !******************************************************************************* 284 WRITE(lunout,*) 'avant q_sat' 285 CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat) 286 WRITE(lunout,*) 'apres q_sat' 287 WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:)) 288 ! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) 289 290 x = 'q'; qd (:,:,:) = 0.0 291 CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,qsat,qd,0.0,jjm, rlonu,rlatv,ib) 292 q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:) 293 294 !--- OZONE CLIMATOLOGY 295 IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz) 296 297 x = 'tsol'; tsol(:) = 0.0 298 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,tsol, 0.0,jjm,rlonu,rlatv,ib) 299 300 x = 'qsol'; qsol(:) = 0.0 301 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,qsol, 0.0,jjm,rlonu,rlatv,ib) 302 303 x = 'snow'; sn(:) = 0.0 304 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,sn, 0.0,jjm,rlonu,rlatv,ib) 305 306 x = 'rads'; radsol(:) = 0.0 307 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib) 308 309 x = 'rugmer'; rugmer(:) = 0.0 310 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib) 311 312 x = 'zmea'; zmea(:) = 0.0 313 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zmea, 0.0,jjm,rlonu,rlatv,ib) 314 315 x = 'zstd'; zstd(:) = 0.0 316 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zstd, 0.0,jjm,rlonu,rlatv,ib) 317 318 x = 'zsig'; zsig(:) = 0.0 319 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zsig, 0.0,jjm,rlonu,rlatv,ib) 320 321 x = 'zgam'; zgam(:) = 0.0 322 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zgam, 0.0,jjm,rlonu,rlatv,ib) 323 324 x = 'zthe'; zthe(:) = 0.0 325 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zthe, 0.0,jjm,rlonu,rlatv,ib) 326 327 x = 'zpic'; zpic(:) = 0.0 328 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zpic, 0.0,jjm,rlonu,rlatv,ib) 329 330 x = 'zval'; zval(:) = 0.0 331 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zval, 0.0,jjm,rlonu,rlatv,ib) 332 333 ! WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273) 334 335 ! Soil ice file reading for soil fraction and soil ice fraction 336 !******************************************************************************* 337 CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid) 338 ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic)) 339 ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic)) 340 ALLOCATE( fraclic(iml_lic,jml_lic)) 341 CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, & 342 lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) 343 CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic) 344 CALL flinclo(fid) 345 346 ! Interpolation on model T-grid 347 !******************************************************************************* 348 WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic 349 ! conversion if coordinates are in degrees 350 IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. 351 IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180. 352 dlon_lic(:)=lon_lic(:,1) 353 dlat_lic(:)=lat_lic(1,:) 354 CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1, & 355 rlonv, rlatu, flic_tmp(1:iim,:) ) 356 flic_tmp(iip1,:)=flic_tmp(1,:) 357 358 !--- To the physical grid 359 CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic)) 360 361 !--- Adequation with soil/sea mask 362 WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 363 WHERE(zmasq(:)<EPSFRA) pctsrf(:,is_lic)=0. 364 pctsrf(:,is_ter)=zmasq(:) 365 DO ji=1,klon 366 IF(zmasq(ji)>EPSFRA) THEN 367 IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN 368 pctsrf(ji,is_lic)=zmasq(ji) 369 pctsrf(ji,is_ter)=0. 370 ELSE 371 pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic) 372 IF(pctsrf(ji,is_ter)<EPSFRA) THEN 373 pctsrf(ji,is_ter)=0. 374 pctsrf(ji,is_lic)=zmasq(ji) 375 END IF 376 END IF 377 END IF 378 END DO 379 380 ! sub-surface ocean and sea ice (sea ice set to zero for start) 381 !******************************************************************************* 382 pctsrf(:,is_oce)=(1.-zmasq(:)) 383 WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0. 384 IF(couple) pctsrf(:,is_oce)=ocemask_fi(:) 385 isst=0 386 WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1 387 388 ! It is checked that the sub-surfaces sum is equal to 1 389 !******************************************************************************* 390 ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA) 391 IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points' 392 CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d) 393 ! WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt) 394 ! WRITE(lunout,*)'zmasq = ' 395 ! WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d) 396 CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 397 WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt) 398 WRITE(lunout,*) 'MASQUE construit : Masque' 399 WRITE(lunout,TRIM(fmt))NINT(masque(:,:)) 400 401 ! Intermediate computation 402 !******************************************************************************* 403 CALL massdair(p3d,masse) 404 WRITE(lunout,*)' ALPHAX ',alphax 405 DO l=1,llm 406 xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l) 407 xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l) 408 xpn=SUM(xppn)/apoln 409 xps=SUM(xpps)/apols 410 masse(:,1 ,l)=xpn 411 masse(:,jjp1,l)=xps 412 END DO 413 q3d(iip1,:,:,:)=q3d(1,:,:,:) 414 phis(iip1,:) = phis(1,:) 415 416 IF(.NOT.letat0) RETURN 417 418 ! Writing 419 !******************************************************************************* 420 CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp) 421 WRITE(lunout,*)'sortie inidissip' 422 itau=0 423 itau_dyn=0 424 itau_phy=0 425 iday=dayref+itau/day_step 426 time=FLOAT(itau-(iday-dayref)*day_step)/day_step 427 IF(time>1.) THEN 428 time=time-1 429 iday=iday+1 430 END IF 431 day_ref=dayref 432 annee_ref=anneeref 433 434 CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) 435 WRITE(lunout,*)'sortie geopot' 436 437 CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & 438 phi, w, pbaru, pbarv, time+iday-dayref) 439 WRITE(lunout,*)'sortie caldyn0' 440 CALL dynredem0( "start.nc", dayref, phis) 441 WRITE(lunout,*)'sortie dynredem0' 442 CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 443 WRITE(lunout,*)'sortie dynredem1' 444 445 ! Physical initial state writting 446 !******************************************************************************* 447 WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad 448 phystep = dtvr * FLOAT(iphysiq) 449 radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) 450 WRITE(lunout,*)'phystep =', phystep, radpas 451 452 ! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs 453 !******************************************************************************* 454 DO i=1,nbsrf; ftsol(:,i) = tsol; END DO 455 DO i=1,nbsrf; snsrf(:,i) = sn; END DO 456 falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6 457 falb1(:,is_oce) = 0.5; falb1(:,is_sic) = 0.6 458 falb2 = falb1 459 evap(:,:) = 0. 460 DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO 461 DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO 462 rain_fall = 0.; snow_fall = 0. 463 solsw = 165.; sollw = -53. 464 t_ancien = 273.15 465 q_ancien = 0. 466 agesno = 0. 467 frugs(:,is_oce) = rugmer(:) 468 frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) 469 frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) 470 frugs(:,is_sic) = 0.001 471 fder = 0.0 472 clwcon = 0.0 473 rnebcon = 0.0 474 ratqs = 0.0 475 run_off_lic_0 = 0.0 476 rugoro = 0.0 477 478 ! Before phyredem calling, surface modules and values to be saved in startphy.nc 479 ! are initialized 480 !******************************************************************************* 481 dummy = 1.0 482 pbl_tke(:,:,:) = 1.e-8 483 zmax0(:) = 40. 484 f0(:) = 1.e-5 485 ema_work1(:,:) = 0. 486 ema_work2(:,:) = 0. 487 wake_deltat(:,:) = 0. 488 wake_deltaq(:,:) = 0. 489 wake_s(:) = 0. 490 wake_cstar(:) = 0. 491 wake_fip(:) = 0. 492 CALL fonte_neige_init(run_off_lic_0) 493 CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil ) 494 CALL phyredem( "startphy.nc" ) 495 496 WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' 497 WRITE(lunout,*)'entree histclo' 498 CALL histclo() 786 499 787 500 #endif 788 501 !#endif of #ifdef CPP_EARTH 789 RETURN 790 ! 791 END SUBROUTINE etat0_netcdf 502 RETURN 503 504 END SUBROUTINE etat0_netcdf 505 ! 506 !------------------------------------------------------------------------------- -
LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90
r1318 r1319 1 1 ! 2 2 ! $Id$ 3 ! 4 C 5 C 6 SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque) 3 !------------------------------------------------------------------------------- 4 ! 5 SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque) 6 ! 7 !------------------------------------------------------------------------------- 8 ! Author : L. Fairhead, 27/01/94 9 !------------------------------------------------------------------------------- 10 ! Purpose: Boundary conditions files building for new model using climatologies. 11 ! Both grids have to be regular. 12 !------------------------------------------------------------------------------- 13 ! Note: This routine is designed to work for Earth 14 !------------------------------------------------------------------------------- 15 ! Modification history: 16 ! * 23/03/1994: Z. X. Li 17 ! * 09/1999: L. Fairhead (netcdf reading in LMDZ.3.3) 18 ! * 07/2001: P. Le Van 19 ! * 11/2009: L. Guez (ozone day & night climatos, see etat0_netcdf.F90) 20 ! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs) 21 !------------------------------------------------------------------------------- 7 22 #ifdef CPP_EARTH 8 ! This routine is designed to work for Earth 9 USE dimphy 10 use phys_state_var_mod , ONLY : pctsrf 11 IMPLICIT none 12 c 13 c------------------------------------------------------------- 14 C Author : L. Fairhead 15 C Date : 27/01/94 16 C Objet : Construction des fichiers de conditions aux limites 17 C pour le nouveau 18 C modele a partir de fichiers de climatologie. Les deux 19 C grilles doivent etre regulieres 20 c 21 c Modifie par z.x.li (le23mars1994) 22 c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999 23 c pour lecture netcdf dans LMDZ.3.3 24 c Modifie par P;Le Van , juillet 2001 25 c------------------------------------------------------------- 26 c 23 USE dimphy 24 USE ioipsl, ONLY : ioget_year_len 25 USE phys_state_var_mod, ONLY : pctsrf 26 USE netcdf, ONLY : NF90_OPEN, NF90_CREATE, NF90_CLOSE, & 27 NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT, & 28 NF90_NOERR, NF90_NOWRITE, NF90_DOUBLE, NF90_GLOBAL, & 29 NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED 30 #endif 31 IMPLICIT NONE 32 !------------------------------------------------------------------------------- 33 ! Arguments: 27 34 #include "dimensions.h" 28 35 #include "paramet.h" 36 #include "iniprint.h" 37 LOGICAL, INTENT(IN) :: interbar ! barycentric interpolation 38 LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag 39 LOGICAL, INTENT(IN) :: oldice ! old way ice computation 40 REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask 41 #ifndef CPP_EARTH 42 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' 43 #else 44 !------------------------------------------------------------------------------- 45 ! Local variables: 29 46 #include "control.h" 30 47 #include "logic.h" 31 #include "netcdf.inc"32 48 #include "comvert.h" 33 49 #include "comgeom2.h" 34 50 #include "comconst.h" 35 cy#include "dimphy.h" 51 #include "indicesol.h" 52 53 !--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) -------- 54 LOGICAL, PARAMETER :: fracterre=.TRUE. 55 56 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 57 CHARACTER(LEN=25) :: icefile, sstfile, dumstr 58 CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', & 59 famipsic='amipbc_sic_1x1.nc ', & 60 fclimsst='amipbc_sst_1x1_clim.nc ', & 61 fclimsic='amipbc_sic_1x1_clim.nc ', & 62 fcpldsst='cpl_atm_sst.nc ', & 63 fcpldsic='cpl_atm_sic.nc ', & 64 frugo ='Rugos.nc ', & 65 falbe ='Albedo.nc ' 66 67 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ 68 REAL, DIMENSION(klon) :: fi_ice, verif 69 REAL, DIMENSION(:,:), POINTER :: phy_rug=>NULL(), phy_ice=>NULL() 70 REAL, DIMENSION(:,:), POINTER :: phy_sst=>NULL(), phy_alb=>NULL() 71 REAL, DIMENSION(:,:), ALLOCATABLE :: phy_bil 72 REAL, DIMENSION(:,:,:), ALLOCATABLE :: pctsrf_t 73 INTEGER :: nbad 74 75 !--- VARIABLES FOR OUTPUT FILE WRITING ----------------------------------------- 76 INTEGER :: ierr, nid, ndim, ntim, k 77 INTEGER, DIMENSION(2) :: dims 78 INTEGER :: id_tim, id_SST, id_BILS, id_RUG, id_ALB 79 INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC 80 INTEGER :: NF90_FORMAT 81 LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files 82 83 !--- MICS (PHYSICAL KEYS, TIME) ------------------------------------------------ 84 ! INTEGER, PARAMETER :: longcles=20 85 ! REAL, DIMENSION(longcles) :: clesphy0 86 INTEGER :: ndays 87 88 !--- INITIALIZATIONS ----------------------------------------------------------- 89 #ifdef NC_DOUBLE 90 NF90_FORMAT=NF90_DOUBLE 91 #else 92 NF90_FORMAT=NF90_FLOAT 93 #endif 94 95 ! CALL conf_gcm(99, .TRUE., clesphy0) 96 pi = 4.*ATAN(1.) 97 rad = 6371229. 98 daysec= 86400. 99 omeg = 2.*pi/daysec 100 g = 9.8 101 kappa = 0.2857143 102 cpp = 1004.70885 103 dtvr = daysec/FLOAT(day_step) 104 CALL inigeom 105 106 !--- Beware: anneeref (from gcm.def) is used to determine output time sampling 107 ndays=ioget_year_len(anneeref) 108 109 !--- RUGOSITY TREATMENT -------------------------------------------------------- 110 WRITE(lunout,*) 'Traitement de la rugosite' 111 CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 112 113 !--- OCEAN TREATMENT ----------------------------------------------------------- 114 PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE. 115 116 ! Input SIC file selection 117 icefile='fake' 118 IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic) 119 IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic) 120 IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic) 121 SELECT CASE(icefile) 122 CASE(famipsic); dumstr='Amip.' 123 CASE(fclimsic); dumstr='Amip climatologique.' 124 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 125 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1) 126 END SELECT 127 ierr=NF90_CLOSE(nid) 128 WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr) 129 130 CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice) 131 132 ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) 133 DO k=1,ndays 134 fi_ice=phy_ice(:,k) 135 WHERE(fi_ice>=1.0 ) fi_ice=1.0 136 WHERE(fi_ice<EPSFRA) fi_ice=0.0 137 IF(fracterre) THEN 138 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 139 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 140 IF(lCPL) THEN ! SIC=pICE*(1-LIC-TER) 141 pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 142 ELSE ! SIC=pICE-LIC 143 pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k) 144 END IF 145 WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0. 146 WHERE(1.0-zmasq<EPSFRA) 147 pctsrf_t(:,is_sic,k)=0.0 148 pctsrf_t(:,is_oce,k)=0.0 149 ELSEWHERE 150 WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq) 151 pctsrf_t(:,is_sic,k)=1.0-zmasq 152 pctsrf_t(:,is_oce,k)=0.0 153 ELSEWHERE 154 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) 155 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA) 156 pctsrf_t(:,is_oce,k)=0.0 157 pctsrf_t(:,is_sic,k)=1.0-zmasq 158 END WHERE 159 END WHERE 160 END WHERE 161 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 162 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 163 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 164 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 165 ELSE 166 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) 167 WHERE(NINT(pctsrf(:,is_ter))==1) 168 pctsrf_t(:,is_sic,k)=0. 169 pctsrf_t(:,is_oce,k)=0. 170 WHERE(fi_ice>=1.e-5) 171 pctsrf_t(:,is_lic,k)=fi_ice 172 pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k) 173 ELSEWHERE 174 pctsrf_t(:,is_lic,k)=0.0 175 END WHERE 176 ELSEWHERE 177 pctsrf_t(:,is_lic,k) = 0.0 178 WHERE(fi_ice>=1.e-5) 179 pctsrf_t(:,is_ter,k)=0.0 180 pctsrf_t(:,is_sic,k)=fi_ice 181 pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k) 182 ELSEWHERE 183 pctsrf_t(:,is_sic,k)=0.0 184 pctsrf_t(:,is_oce,k)=1.0 185 END WHERE 186 END WHERE 187 verif=sum(pctsrf_t(:,:,k),dim=2) 188 nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5) 189 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 190 END IF 191 END DO 192 DEALLOCATE(phy_ice) 193 194 !--- SST TREATMENT ------------------------------------------------------------- 195 WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE. 196 197 ! Input SST file selection 198 sstfile='fake' 199 IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst) 200 IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst) 201 IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst) 202 SELECT CASE(icefile) 203 CASE(famipsic); dumstr='Amip.' 204 CASE(fclimsic); dumstr='Amip climatologique.' 205 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 206 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1) 207 END SELECT 208 ierr=NF90_CLOSE(nid) 209 WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr) 210 211 CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap) 212 213 !--- ALBEDO TREATMENT ---------------------------------------------------------- 214 WRITE(lunout,*) 'Traitement de l albedo' 215 CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb) 216 217 !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- 218 ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0 219 220 !--- OUTPUT FILE WRITING ------------------------------------------------------- 221 WRITE(lunout,*) 'Ecriture du fichier limit' 222 223 !--- File creation 224 ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid) 225 ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites") 226 227 !--- Dimensions creation 228 ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim) 229 ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim) 230 231 dims=(/ndim,ntim/) 232 233 !--- Variables creation 234 ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,dims,id_tim) 235 ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE) 236 ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC) 237 ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER) 238 ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC) 239 ierr=NF90_DEF_VAR(nid,"SST", NF90_FORMAT,dims,id_SST) 240 ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS) 241 ierr=NF90_DEF_VAR(nid,"ALB", NF90_FORMAT,dims,id_ALB) 242 ierr=NF90_DEF_VAR(nid,"RUG", NF90_FORMAT,dims,id_RUG) 243 244 !--- Attributes creation 245 ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee") 246 ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean") 247 ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer") 248 ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre") 249 ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice") 250 ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer") 251 ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol") 252 ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface") 253 ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite") 254 255 ierr=NF90_ENDDEF(nid) 256 257 !--- Variables saving 258 ierr=NF90_PUT_VAR(nid,id_tim,(/(DBLE(k),k=1,ndays)/)) 259 ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/)) 260 ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/)) 261 ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/)) 262 ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/)) 263 ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/)) 264 ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/)) 265 ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/)) 266 ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/)) 267 268 ierr=NF90_CLOSE(nid) 269 270 DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) 271 #endif 272 ! of #ifdef CPP_EARTH 273 STOP 274 275 276 !=============================================================================== 277 ! 278 CONTAINS 279 ! 280 !=============================================================================== 281 282 283 !------------------------------------------------------------------------------- 284 ! 285 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask) 286 ! 287 !------------------------------------------------------------------------------- 288 ! Comments: 289 ! There are two assumptions concerning the NetCDF files, that are satisfied 290 ! with files that are conforming NC convention: 291 ! 1) The last dimension of the variables used is the time record. 292 ! 2) Dimensional variables have the same names as corresponding dimensions. 293 !------------------------------------------------------------------------------- 294 USE netcdf, ONLY : NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, & 295 NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, & 296 NF90_GET_VAR, NF90_GET_ATT 297 USE dimphy, ONLY : klon 298 USE phys_state_var_mod, ONLY : pctsrf 299 IMPLICIT NONE 300 #include "dimensions.h" 301 #include "paramet.h" 302 #include "comgeom2.h" 303 #include "control.h" 36 304 #include "indicesol.h" 37 305 #include "iniprint.h" 38 c 39 c----------------------------------------------------------------------- 40 LOGICAL interbar, extrap, oldice 41 42 REAL phy_nat(klon,360), phy_nat0(klon) 43 REAL phy_alb(klon,360) 44 REAL phy_sst(klon,360) 45 REAL phy_bil(klon,360) 46 REAL phy_rug(klon,360) 47 REAL phy_ice(klon) 48 c 49 real pctsrf_t(klon,nbsrf,360) 50 51 REAL verif 52 53 REAL masque(iip1,jjp1) 54 REAL mask(iim,jjp1) 55 CPB 56 C newlmt indique l'utilisation de la sous-maille fractionnelle 57 C tandis que l'ancien codage utilise l'indicateur du sol (0,1,2,3) 58 LOGICAL newlmt, fracterre 59 PARAMETER(newlmt=.TRUE.) 60 PARAMETER(fracterre = .TRUE.) 61 62 C Declarations pour le champ de depart 63 INTEGER imdep, jmdep,lmdep 64 INTEGER tbid 65 PARAMETER ( tbid = 60 ) ! >52 semaines 66 REAL timecoord(tbid) 67 c 68 REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:) 69 REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:) 70 REAL , ALLOCATABLE :: dlon(:), dlat(:) 71 REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:) 72 REAL , ALLOCATABLE :: champ_msk(:), champ(:) 73 REAL , ALLOCATABLE :: work(:,:) 74 75 CHARACTER*25 title 76 77 C Declarations pour le champ interpole 2D 78 REAL champint(iim,jjp1) 79 real chmin,chmax 80 81 C Declarations pour le champ interpole 3D 82 REAL champtime(iim,jjp1,tbid) 83 REAL timeyear(tbid) 84 REAL champan(iip1,jjp1,366) 85 86 C Declarations pour l'inteprolation verticale 87 REAL ax(tbid), ay(tbid) 88 REAL by 89 REAL yder(tbid) 90 91 92 INTEGER ierr 93 INTEGER dimfirst(3) 94 INTEGER dimlast(3) 95 c 96 INTEGER nid, ndim, ntim 97 INTEGER dims(2), debut(2), epais(2) 98 INTEGER id_tim 99 INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB 100 CPB 101 INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC 102 103 INTEGER i, j, k, l, ji 104 c declarations pour lecture glace de mer 105 INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret 106 INTEGER :: itaul(1), fid 107 REAL :: lev(1), date, dt 108 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic 109 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic 110 REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic 111 REAL :: flic_tmp(iip1, jjp1) 112 113 c Diverses variables locales 114 REAL time 115 ! pour la lecture du fichier masque ocean 116 integer :: nid_o2a 117 logical :: couple = .false. 118 INTEGER :: iml_omask, jml_omask 119 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask 120 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask 121 REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp 122 real, dimension(klon) :: ocemask_fi 123 124 INTEGER longcles 125 PARAMETER ( longcles = 20 ) 126 REAL clesphy0 ( longcles ) 127 #include "serre.h" 128 INTEGER ncid,varid,ndimid(4),dimid 129 character*30 namedim 130 CHARACTER*80 :: varname 131 132 cIM28/02/2002 <== PM 133 REAL tmidmonth(12) 134 SAVE tmidmonth 135 DATA tmidmonth/15,45,75,105,135,165,195,225,255,285,315,345/ 136 137 c initialisations: 138 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 139 140 141 pi = 4. * ATAN(1.) 142 rad = 6 371 229. 143 omeg = 4.* ASIN(1.)/(24.*3600.) 144 g = 9.8 145 daysec = 86400. 146 kappa = 0.2857143 147 cpp = 1004.70885 148 dtvr = daysec/FLOAT(day_step) 149 CALL inigeom 150 c 151 C Traitement du relief au sol 152 c 153 write(*,*) 'Traitement du relief au sol pour fabriquer masque' 154 ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid) 155 156 if (ierr.ne.0) then 157 print *, NF_STRERROR(ierr) 158 STOP 159 ENDIF 160 161 ierr = NF_INQ_VARID(ncid,'RELIEF',varid) 162 if (ierr.ne.0) then 163 print *, NF_STRERROR(ierr) 164 STOP 165 ENDIF 166 ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) 167 if (ierr.ne.0) then 168 print *, NF_STRERROR(ierr) 169 STOP 170 ENDIF 171 ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) 172 if (ierr.ne.0) then 173 print *, NF_STRERROR(ierr) 174 STOP 175 ENDIF 176 print*,'variable ', namedim, 'dimension ', imdep 177 ierr = NF_INQ_VARID(ncid,namedim,dimid) 178 if (ierr.ne.0) then 179 print *, NF_STRERROR(ierr) 180 STOP 181 ENDIF 182 183 ALLOCATE( lonmsk_ini(imdep) ) 184 ALLOCATE( dlon_msk(imdep) ) 185 186 #ifdef NC_DOUBLE 187 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini) 188 #else 189 ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini) 190 #endif 191 192 c 193 if (ierr.ne.0) then 194 print *, NF_STRERROR(ierr) 195 STOP 196 ENDIF 197 ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) 198 if (ierr.ne.0) then 199 print *, NF_STRERROR(ierr) 200 STOP 201 ENDIF 202 print*,'variable ', namedim, 'dimension ', jmdep 203 ierr = NF_INQ_VARID(ncid,namedim,dimid) 204 if (ierr.ne.0) then 205 print *, NF_STRERROR(ierr) 206 STOP 207 ENDIF 208 209 ALLOCATE( latmsk_ini(jmdep) ) 210 ALLOCATE( dlat_msk(jmdep) ) 211 ALLOCATE( champ_msk(imdep*jmdep) ) 212 213 #ifdef NC_DOUBLE 214 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini) 215 #else 216 ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini) 217 #endif 218 c 219 if (ierr.ne.0) then 220 print *, NF_STRERROR(ierr) 221 STOP 222 ENDIF 223 #ifdef NC_DOUBLE 224 ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk) 225 #else 226 ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk) 227 #endif 228 c 229 if (ierr.ne.0) then 230 print *, NF_STRERROR(ierr) 231 STOP 232 ENDIF 233 c 234 title='RELIEF' 235 236 CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini, 237 . dlon_msk, dlat_msk, champ_msk, interbar ) 238 239 DO i = 1, iim 240 DO j = 1, jjp1 241 mask(i,j) = masque(i,j) 242 ENDDO 243 ENDDO 244 WRITE(*,*) 'MASK:' 245 WRITE(*,'(96i1)')INT(mask) 246 ierr = NF_CLOSE(ncid) 247 c 248 c 249 C Traitement de la rugosite 250 c 251 PRINT*, 'Traitement de la rugosite' 252 ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid) 253 if (ierr.ne.0) then 254 print *, NF_STRERROR(ierr) 255 STOP 256 ENDIF 257 258 ierr = NF_INQ_VARID(ncid,'RUGOS',varid) 259 if (ierr.ne.0) then 260 print *, NF_STRERROR(ierr) 261 STOP 262 ENDIF 263 ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) 264 if (ierr.ne.0) then 265 print *, NF_STRERROR(ierr) 266 STOP 267 ENDIF 268 ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) 269 if (ierr.ne.0) then 270 print *, NF_STRERROR(ierr) 271 STOP 272 ENDIF 273 print*,'variable ', namedim, 'dimension ', imdep 274 ierr = NF_INQ_VARID(ncid,namedim,dimid) 275 if (ierr.ne.0) then 276 print *, NF_STRERROR(ierr) 277 STOP 278 ENDIF 279 280 ALLOCATE( dlon_ini(imdep) ) 281 ALLOCATE( dlon(imdep) ) 282 283 #ifdef NC_DOUBLE 284 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) 285 #else 286 ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) 287 #endif 288 if (ierr.ne.0) then 289 print *, NF_STRERROR(ierr) 290 STOP 291 ENDIF 292 ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) 293 if (ierr.ne.0) then 294 print *, NF_STRERROR(ierr) 295 STOP 296 ENDIF 297 print*,'variable ', namedim, 'dimension ', jmdep 298 ierr = NF_INQ_VARID(ncid,namedim,dimid) 299 if (ierr.ne.0) then 300 print *, NF_STRERROR(ierr) 301 STOP 302 ENDIF 303 304 ALLOCATE( dlat_ini(jmdep) ) 305 ALLOCATE( dlat(jmdep) ) 306 307 #ifdef NC_DOUBLE 308 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) 309 #else 310 ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) 311 #endif 312 if (ierr.ne.0) then 313 print *, NF_STRERROR(ierr) 314 STOP 315 ENDIF 316 ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) 317 if (ierr.ne.0) then 318 print *, NF_STRERROR(ierr) 319 STOP 320 ENDIF 321 print*,'variable ', namedim, 'dimension ', lmdep 322 ierr = NF_INQ_VARID(ncid,namedim,dimid) 323 if (ierr.ne.0) then 324 print *, NF_STRERROR(ierr) 325 STOP 326 ENDIF 327 #ifdef NC_DOUBLE 328 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) 329 #else 330 ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) 331 #endif 332 if (ierr.ne.0) then 333 print *, NF_STRERROR(ierr) 334 STOP 335 ENDIF 336 c 337 ALLOCATE( champ(imdep*jmdep) ) 338 339 DO 200 l = 1, lmdep 340 dimfirst(1) = 1 341 dimfirst(2) = 1 342 dimfirst(3) = l 343 c 344 dimlast(1) = imdep 345 dimlast(2) = jmdep 346 dimlast(3) = 1 347 c 348 PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) 349 print*,dimfirst,dimlast 350 #ifdef NC_DOUBLE 351 ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) 352 #else 353 ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) 354 #endif 355 if (ierr.ne.0) then 356 print *, NF_STRERROR(ierr) 357 STOP 358 ENDIF 359 360 title = 'Rugosite Amip ' 361 c 362 CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, 363 . dlon, dlat, champ, interbar ) 364 365 IF ( interbar ) THEN 366 DO j = 1, imdep * jmdep 367 champ(j) = LOG(champ(j)) 368 ENDDO 369 370 IF( l.EQ.1 ) THEN 371 WRITE(6,*) '-------------------------------------------------', 372 ,'------------------------' 373 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 374 , ' pour la rugosite $$$ ' 375 WRITE(6,*) '-------------------------------------------------', 376 ,'------------------------' 377 ENDIF 378 CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ , 379 , iim,jjm,rlonu,rlatv, jjp1,champint ) 380 DO j=1,jjp1 381 DO i=1,iim 382 champint(i,j)=EXP(champint(i,j)) 383 ENDDO 384 ENDDO 385 386 DO j = 1, jjp1 387 DO i = 1, iim 388 IF(NINT(mask(i,j)).NE.1) THEN 389 champint( i,j ) = 0.001 390 ENDIF 391 ENDDO 392 ENDDO 393 ELSE 394 CALL rugosite(imdep, jmdep, dlon, dlat, champ, 395 . iim, jjp1, rlonv, rlatu, champint, mask) 396 ENDIF 397 DO j = 1,jjp1 398 DO i = 1, iim 399 champtime (i,j,l) = champint(i,j) 400 ENDDO 401 ENDDO 402 200 CONTINUE 403 c 404 DO l = 1, lmdep 405 timeyear(l) = timecoord(l) 406 ENDDO 407 408 PRINT 222, timeyear(:lmdep) 409 222 FORMAT(2x,' Time year ',10f6.1) 410 c 411 412 PRINT*, 'Interpolation temporelle dans l annee' 413 414 DO j = 1, jjp1 415 DO i = 1, iim 416 DO l = 1, lmdep 417 ax(l) = timeyear(l) 418 ay(l) = champtime (i,j,l) 419 ENDDO 420 CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) 421 DO k = 1, 360 422 time = FLOAT(k-1) 423 CALL SPLINT(ax,ay,yder,lmdep,time,by) 424 champan(i,j,k) = by 425 ENDDO 426 ENDDO 427 ENDDO 428 DO k = 1, 360 429 DO j = 1, jjp1 430 champan(iip1,j,k) = champan(1,j,k) 431 ENDDO 432 IF ( k.EQ.10 ) THEN 433 DO j = 1, jjp1 434 CALL minmax( iip1,champan(1,j,10),chmin,chmax ) 435 PRINT *,' Rugosite au temps 10 ', chmin,chmax,j 436 ENDDO 437 ENDIF 438 ENDDO 439 c 440 DO k = 1, 360 441 CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k)) 442 ENDDO 443 c 444 ierr = NF_CLOSE(ncid) 445 446 DEALLOCATE( dlon ) 447 DEALLOCATE( dlon_ini ) 448 DEALLOCATE( dlat ) 449 DEALLOCATE( dlat_ini ) 450 DEALLOCATE( champ ) 451 c 452 c 453 C Traitement de la glace oceanique 454 c 455 PRINT*, 'Traitement de la glace oceanique' 456 457 ierr = NF_OPEN('amipbc_sic_1x1.nc', NF_NOWRITE, ncid) 458 if (ierr.ne.0) THEN 459 ierr = NF_OPEN('amipbc_sic_1x1_clim.nc', NF_NOWRITE, ncid) 460 if (ierr.ne.0) THEN 461 print *, NF_STRERROR(ierr) 462 STOP 463 endif 464 ENDIF 465 466 cIM22/02/2002 467 cIM07/03/2002 AMIP.nc & amip79to95.nc 468 cIM ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid) 469 cIM07/03/2002 amipbc_sic_1x1_clim.nc & amipbc_sic_1x1.nc 470 ierr = NF_INQ_VARID(ncid,'sicbcs',varid) 471 cIM22/02/2002 472 if (ierr.ne.0) then 473 print *, NF_STRERROR(ierr),'sicbcs' 474 STOP 475 ENDIF 476 ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) 477 if (ierr.ne.0) then 478 print *, NF_STRERROR(ierr) 479 STOP 480 ENDIF 481 ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) 482 if (ierr.ne.0) then 483 print *, NF_STRERROR(ierr) 484 STOP 485 ENDIF 486 print*,'variable ', namedim, 'dimension ', imdep 487 ierr = NF_INQ_VARID(ncid,namedim,dimid) 488 if (ierr.ne.0) then 489 print *, NF_STRERROR(ierr) 490 STOP 491 ENDIF 492 493 ALLOCATE ( dlon_ini(imdep) ) 494 ALLOCATE ( dlon(imdep) ) 495 496 #ifdef NC_DOUBLE 497 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) 498 #else 499 ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) 500 #endif 501 if (ierr.ne.0) then 502 print *, NF_STRERROR(ierr) 503 STOP 504 ENDIF 505 ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) 506 if (ierr.ne.0) then 507 print *, NF_STRERROR(ierr) 508 STOP 509 ENDIF 510 print*,'variable ', namedim, jmdep 511 ierr = NF_INQ_VARID(ncid,namedim,dimid) 512 if (ierr.ne.0) then 513 print *, NF_STRERROR(ierr) 514 STOP 515 ENDIF 516 517 ALLOCATE ( dlat_ini(jmdep) ) 518 ALLOCATE ( dlat(jmdep) ) 519 520 #ifdef NC_DOUBLE 521 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) 522 #else 523 ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) 524 #endif 525 if (ierr.ne.0) then 526 print *, NF_STRERROR(ierr) 527 STOP 528 ENDIF 529 ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) 530 if (ierr.ne.0) then 531 print *, NF_STRERROR(ierr) 532 STOP 533 ENDIF 534 print*,'variable ', namedim, lmdep 535 cIM28/02/2002 536 cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours 537 c Ici on suppose qu'on a 12 mois (de 30 jours). 538 IF (lmdep.NE.12) THEN 539 print *, 'Unknown AMIP file: not 12 months ?' 540 STOP 541 ENDIF 542 cIM28/02/2002 543 ierr = NF_INQ_VARID(ncid,namedim,dimid) 544 if (ierr.ne.0) then 545 print *, NF_STRERROR(ierr) 546 STOP 547 ENDIF 548 #ifdef NC_DOUBLE 549 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) 550 #else 551 ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) 552 #endif 553 if (ierr.ne.0) then 554 print *, NF_STRERROR(ierr) 555 STOP 556 ENDIF 557 c 558 ALLOCATE ( champ(imdep*jmdep) ) 559 560 DO l = 1, lmdep 561 dimfirst(1) = 1 562 dimfirst(2) = 1 563 dimfirst(3) = l 564 c 565 dimlast(1) = imdep 566 dimlast(2) = jmdep 567 dimlast(3) = 1 568 c 569 PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) 570 #ifdef NC_DOUBLE 571 ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) 572 #else 573 ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) 574 #endif 575 if (ierr.ne.0) then 576 print *, NF_STRERROR(ierr) 577 STOP 578 ENDIF 579 580 title = 'Sea-ice Amip ' 581 c 582 CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, 583 . dlon, dlat, champ, interbar ) 584 c 585 IF( oldice ) THEN 586 CALL sea_ice(imdep, jmdep, dlon, dlat, champ, 587 . iim, jjp1, rlonv, rlatu, champint ) 588 ELSEIF ( interbar ) THEN 589 IF( l.EQ.1 ) THEN 590 WRITE(6,*) '-------------------------------------------------', 591 ,'------------------------' 592 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 593 , ' pour Sea-ice Amip $$$ ' 594 WRITE(6,*) '-------------------------------------------------', 595 ,'------------------------' 596 ENDIF 597 cIM07/03/2002 598 cIM22/02/2002 : Sea-ice Amip entre 0. et 1. 599 cIM PRINT*,'SUB. limit_netcdf.F IM : Sea-ice et SST Amip_new clim' 600 cIM DO j = 1, imdep * jmdep 601 cIM28/02/2002 <==PM champ(j) = champ(j)/100. 602 cIM14/03/2002 champ(j) = max(0.0,(min(1.0, (champ(j)/100.) ))) 603 cIM champ(j) = amax1(0.0,(amin1(1.0, (champ(j)/100.) ))) 604 cIM ENDDO 605 cIM22/02/2002 606 CALL inter_barxy ( imdep,jmdep -1,dlon, dlat , 607 , champ, iim, jjm, rlonu, rlatv, jjp1, champint ) 608 ELSE 609 CALL sea_ice(imdep, jmdep, dlon, dlat, champ, 610 . iim, jjp1, rlonv, rlatu, champint ) 611 ENDIF 612 DO j = 1,jjp1 613 DO i = 1, iim 614 champtime (i,j,l) = champint(i,j) 615 ENDDO 616 ENDDO 617 ENDDO 618 c 619 DO l = 1, lmdep 620 cIM28/02/2002 <== PM timeyear(l) = timecoord(l) 621 cIM timeyear(l) = timecoord(l) 622 cIM07/03/2002 623 timeyear(l) = tmidmonth(l) 624 ENDDO 625 PRINT 222, timeyear(:lmdep) 626 c 627 PRINT*, 'Interpolation temporelle' 628 DO j = 1, jjp1 629 DO i = 1, iim 630 DO l = 1, lmdep 631 ax(l) = timeyear(l) 632 ay(l) = champtime (i,j,l) 633 ENDDO 634 CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) 635 DO k = 1, 360 636 time = FLOAT(k-1) 637 CALL SPLINT(ax,ay,yder,lmdep,time,by) 638 champan(i,j,k) = by 639 ENDDO 640 ENDDO 641 ENDDO 642 DO k = 1, 360 643 DO j = 1, jjp1 644 champan(iip1, j, k) = champan(1, j, k) 645 ENDDO 646 IF ( k.EQ.10 ) THEN 647 DO j = 1, jjp1 648 CALL minmax( iip1,champan(1,j,10),chmin,chmax ) 649 PRINT *,' Sea ice au temps 10 ', chmin,chmax,j 650 ENDDO 651 ENDIF 652 ENDDO 653 c 654 cIM14/03/2002 : Sea-ice Amip entre 0. et 1. 655 PRINT*,'SUB. limit_netcdf.F IM : Sea-ice Amipbc ' 656 DO k = 1, 360 657 DO j = 1, jjp1 658 DO i = 1, iim 659 champan(i, j, k) = 660 $ amax1(0.0,(amin1(1.0,(champan(i, j, k)/100.)))) 661 ENDDO 662 champan(iip1, j, k) = champan(1, j, k) 663 ENDDO 664 ENDDO 665 cIM14/03/2002 666 667 DO k = 1, 360 668 CALL gr_dyn_fi(1, iip1, jjp1, klon, 669 . champan(1,1,k), phy_ice) 670 IF ( newlmt) THEN 671 672 CPB en attendant de mettre fraction de terre 673 c 674 WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1. 675 WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0. 676 c 677 IF (fracterre ) THEN 678 c WRITE(*,*) 'passe dans cas fracterre' 679 pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter) 680 pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic) 681 pctsrf_t(1:klon,is_sic,k) = phy_ice(1:klon) 682 $ - pctsrf_t(1:klon,is_lic,k) 683 c Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP 684 WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0) 685 pctsrf_t(1:klon,is_sic,k) = 0. 686 END WHERE 687 WHERE( 1. - zmasq(1:klon) .LT. EPSFRA) 688 pctsrf_t(1:klon,is_sic,k) = 0. 689 pctsrf_t(1:klon,is_oce,k) = 0. 690 END WHERE 691 DO i = 1, klon 692 IF ( 1. - zmasq(i) .GT. EPSFRA) THEN 693 IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN 694 pctsrf_t(i,is_sic,k) = 1 - zmasq(i) 695 pctsrf_t(i,is_oce,k) = 0. 696 ELSE 697 pctsrf_t(i,is_oce,k) = 1 - zmasq(i) 698 $ - pctsrf_t(i,is_sic,k) 699 IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN 700 pctsrf_t(i,is_oce,k) = 0. 701 pctsrf_t(i,is_sic,k) = 1 - zmasq(i) 702 ENDIF 703 ENDIF 704 ENDIF 705 if (pctsrf_t(i,is_oce,k) .lt. 0.) then 706 WRITE(*,*) 'pb sous maille au point : i,k ' 707 $ , i,k,pctsrf_t(:,is_oce,k) 708 ENDIF 709 IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) + 710 $ pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k) - 1.) 711 $ .GT. EPSFRA) THEN 712 WRITE(*,*) 'physiq : pb sous surface au point ', i, 713 $ pctsrf_t(i, 1 : nbsrf,k), phy_ice(i) 714 ENDIF 715 END DO 716 ELSE 717 DO i = 1, klon 718 pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter) 719 IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN 720 pctsrf_t(i,is_sic,k) = 0. 721 pctsrf_t(i,is_oce,k) = 0. 722 IF(phy_ice(i) .GE. 1.e-5) THEN 723 pctsrf_t(i,is_lic,k) = phy_ice(i) 724 pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k) 725 . - pctsrf_t(i,is_lic,k) 726 ELSE 727 pctsrf_t(i,is_lic,k) = 0. 728 ENDIF 729 ELSE 730 pctsrf_t(i,is_lic,k) = 0. 731 IF(phy_ice(i) .GE. 1.e-5) THEN 732 pctsrf_t(i,is_ter,k) = 0. 733 pctsrf_t(i,is_sic,k) = phy_ice(i) 734 pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k) 735 ELSE 736 pctsrf_t(i,is_sic,k) = 0. 737 pctsrf_t(i,is_oce,k) = 1. 738 ENDIF 739 ENDIF 740 verif = pctsrf_t(i,is_ter,k) + 741 . pctsrf_t(i,is_oce,k) + 742 . pctsrf_t(i,is_sic,k) + 743 . pctsrf_t(i,is_lic,k) 744 IF ( verif .LT. 1. - 1.e-5 .OR. 745 $ verif .GT. 1 + 1.e-5) THEN 746 WRITE(*,*) 'pb sous maille au point : i,k,verif ' 747 $ , i,k,verif 748 ENDIF 749 END DO 750 ENDIF 751 ELSE 752 DO i = 1, klon 753 phy_nat(i,k) = phy_nat0(i) 754 IF ( (phy_ice(i) - 0.5).GE.1.e-5 ) THEN 755 IF (NINT(phy_nat0(i)).EQ.0) THEN 756 phy_nat(i,k) = 3.0 757 ELSE 758 phy_nat(i,k) = 2.0 759 ENDIF 760 ENDIF 761 IF( NINT(phy_nat(i,k)).EQ.0 ) THEN 762 IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001 763 ENDIF 764 END DO 765 ENDIF 766 ENDDO 767 c 768 769 ierr = NF_CLOSE(ncid) 770 c 771 DEALLOCATE( dlon ) 772 DEALLOCATE( dlon_ini ) 773 DEALLOCATE( dlat ) 774 DEALLOCATE( dlat_ini ) 775 DEALLOCATE( champ ) 776 777 477 continue 778 c 779 C Traitement de la sst 780 c 781 PRINT*, 'Traitement de la sst' 782 c ierr = NF_OPEN('AMIP_SST.nc', NF_NOWRITE, ncid) 783 ierr = NF_OPEN('amipbc_sst_1x1.nc', NF_NOWRITE, ncid) 784 if (ierr.ne.0) THEN 785 ierr = NF_OPEN('amipbc_sst_1x1_clim.nc', NF_NOWRITE, ncid) 786 if (ierr.ne.0) THEN 787 print *, NF_STRERROR(ierr) 788 STOP 789 endif 790 ENDIF 791 792 cIM22/02/2002 793 cIM ierr = NF_INQ_VARID(ncid,'SST',varid) 794 ierr = NF_INQ_VARID(ncid,'tosbcs',varid) 795 cIM22/02/2002 796 if (ierr.ne.0) then 797 print *, NF_STRERROR(ierr) 798 STOP 799 ENDIF 800 ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) 801 if (ierr.ne.0) then 802 print *, NF_STRERROR(ierr) 803 STOP 804 ENDIF 805 ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) 806 if (ierr.ne.0) then 807 print *, NF_STRERROR(ierr) 808 STOP 809 ENDIF 810 print*,'variable SST ', namedim,'dimension ', imdep 811 ierr = NF_INQ_VARID(ncid,namedim,dimid) 812 if (ierr.ne.0) then 813 print *, NF_STRERROR(ierr) 814 STOP 815 ENDIF 816 817 ALLOCATE( dlon_ini(imdep) ) 818 ALLOCATE( dlon(imdep) ) 819 820 #ifdef NC_DOUBLE 821 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) 822 #else 823 ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) 824 #endif 825 826 if (ierr.ne.0) then 827 print *, NF_STRERROR(ierr) 828 STOP 829 ENDIF 830 ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) 831 if (ierr.ne.0) then 832 print *, NF_STRERROR(ierr) 833 STOP 834 ENDIF 835 print*,'variable SST ', namedim, 'dimension ', jmdep 836 ierr = NF_INQ_VARID(ncid,namedim,dimid) 837 if (ierr.ne.0) then 838 print *, NF_STRERROR(ierr) 839 STOP 840 ENDIF 841 842 ALLOCATE( dlat_ini(jmdep) ) 843 ALLOCATE( dlat(jmdep) ) 844 845 #ifdef NC_DOUBLE 846 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) 847 #else 848 ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) 849 #endif 850 if (ierr.ne.0) then 851 print *, NF_STRERROR(ierr) 852 STOP 853 ENDIF 854 ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) 855 if (ierr.ne.0) then 856 print *, NF_STRERROR(ierr) 857 STOP 858 ENDIF 859 print*,'variable ', namedim, 'dimension ', lmdep 860 cIM28/02/2002 861 cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours 862 c Ici on suppose qu'on a 12 mois (de 30 jours). 863 IF (lmdep.NE.12) THEN 864 print *, 'Unknown AMIP file: not 12 months ?' 865 STOP 866 ENDIF 867 cIM28/02/2002 868 ierr = NF_INQ_VARID(ncid,namedim,dimid) 869 if (ierr.ne.0) then 870 print *, NF_STRERROR(ierr) 871 STOP 872 ENDIF 873 #ifdef NC_DOUBLE 874 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) 875 #else 876 ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) 877 #endif 878 if (ierr.ne.0) then 879 print *, NF_STRERROR(ierr) 880 STOP 881 ENDIF 882 883 ALLOCATE( champ(imdep*jmdep) ) 884 IF( extrap ) THEN 885 ALLOCATE ( work(imdep,jmdep) ) 886 ENDIF 887 c 888 DO l = 1, lmdep 889 dimfirst(1) = 1 890 dimfirst(2) = 1 891 dimfirst(3) = l 892 c 893 dimlast(1) = imdep 894 dimlast(2) = jmdep 895 dimlast(3) = 1 896 c 897 PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) 898 #ifdef NC_DOUBLE 899 ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) 900 #else 901 ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) 902 #endif 903 if (ierr.ne.0) then 904 print *, NF_STRERROR(ierr) 905 STOP 906 ENDIF 907 908 title='Sst Amip' 909 c 910 CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, 911 . dlon, dlat, champ, interbar ) 912 IF ( extrap ) THEN 913 CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work) 914 ENDIF 915 c 916 917 IF ( interbar ) THEN 918 IF( l.EQ.1 ) THEN 919 WRITE(6,*) '-------------------------------------------------', 920 ,'------------------------' 921 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 922 , ' pour la Sst Amip $$$ ' 923 WRITE(6,*) '-------------------------------------------------', 924 ,'------------------------' 925 ENDIF 926 CALL inter_barxy ( imdep,jmdep -1,dlon, dlat , 927 , champ, iim, jjm, rlonu, rlatv, jjp1, champint ) 928 ELSE 929 CALL grille_m(imdep, jmdep, dlon, dlat, champ, 930 . iim, jjp1, rlonv, rlatu, champint ) 931 ENDIF 932 933 DO j = 1,jjp1 934 DO i = 1, iim 935 champtime (i,j,l) = champint(i,j) 936 ENDDO 937 ENDDO 938 ENDDO 939 c 940 DO l = 1, lmdep 941 cIM28/02/2002 <==PM timeyear(l) = timecoord(l) 942 timeyear(l) = tmidmonth(l) 943 ENDDO 944 print 222, timeyear(:lmdep) 945 c 946 C interpolation temporelle 947 DO j = 1, jjp1 948 DO i = 1, iim 949 DO l = 1, lmdep 950 ax(l) = timeyear(l) 951 ay(l) = champtime (i,j,l) 952 ENDDO 953 CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) 954 DO k = 1, 360 955 time = FLOAT(k-1) 956 CALL SPLINT(ax,ay,yder,lmdep,time,by) 957 champan(i,j,k) = by 958 ENDDO 959 ENDDO 960 ENDDO 961 DO k = 1, 360 962 DO j = 1, jjp1 963 champan(iip1,j,k) = champan(1,j,k) 964 ENDDO 965 IF ( k.EQ.10 ) THEN 966 DO j = 1, jjp1 967 CALL minmax( iip1,champan(1,j,10),chmin,chmax ) 968 PRINT *,' SST au temps 10 ', chmin,chmax,j 969 ENDDO 970 ENDIF 971 ENDDO 972 c 973 cIM14/03/2002 : SST amipbc greater then 271.38 974 PRINT*,'SUB. limit_netcdf.F IM : SST Amipbc >= 271.38 ' 975 DO k = 1, 360 976 DO j = 1, jjp1 977 DO i = 1, iim 978 champan(i, j, k) = amax1(champan(i, j, k), 271.38) 979 ENDDO 980 champan(iip1, j, k) = champan(1, j, k) 981 ENDDO 982 ENDDO 983 cIM14/03/2002 984 DO k = 1, 360 985 CALL gr_dyn_fi(1, iip1, jjp1, klon, 986 . champan(1,1,k), phy_sst(1,k)) 987 ENDDO 988 c 989 ierr = NF_CLOSE(ncid) 990 c 991 DEALLOCATE( dlon ) 992 DEALLOCATE( dlon_ini ) 993 DEALLOCATE( dlat ) 994 DEALLOCATE( dlat_ini ) 995 DEALLOCATE( champ ) 996 c 997 C Traitement de l'albedo 998 c 999 PRINT*, 'Traitement de l albedo' 1000 ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid) 1001 if (ierr.ne.0) then 1002 print *, NF_STRERROR(ierr) 1003 STOP 1004 ENDIF 1005 ierr = NF_INQ_VARID(ncid,'ALBEDO',varid) 1006 if (ierr.ne.0) then 1007 print *, NF_STRERROR(ierr) 1008 STOP 1009 ENDIF 1010 ierr = NF_INQ_VARDIMID (ncid,varid,ndimid) 1011 if (ierr.ne.0) then 1012 print *, NF_STRERROR(ierr) 1013 STOP 1014 ENDIF 1015 ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep) 1016 if (ierr.ne.0) then 1017 print *, NF_STRERROR(ierr) 1018 STOP 1019 ENDIF 1020 print*,'variable ', namedim, 'dimension ', imdep 1021 ierr = NF_INQ_VARID(ncid,namedim,dimid) 1022 if (ierr.ne.0) then 1023 print *, NF_STRERROR(ierr) 1024 STOP 1025 ENDIF 1026 1027 ALLOCATE ( dlon_ini(imdep) ) 1028 ALLOCATE ( dlon(imdep) ) 1029 1030 #ifdef NC_DOUBLE 1031 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini) 1032 #else 1033 ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini) 1034 #endif 1035 if (ierr.ne.0) then 1036 print *, NF_STRERROR(ierr) 1037 STOP 1038 ENDIF 1039 ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep) 1040 if (ierr.ne.0) then 1041 print *, NF_STRERROR(ierr) 1042 STOP 1043 ENDIF 1044 print*,'variable ', namedim, 'dimension ', jmdep 1045 ierr = NF_INQ_VARID(ncid,namedim,dimid) 1046 if (ierr.ne.0) then 1047 print *, NF_STRERROR(ierr) 1048 STOP 1049 ENDIF 1050 1051 ALLOCATE ( dlat_ini(jmdep) ) 1052 ALLOCATE ( dlat(jmdep) ) 1053 1054 #ifdef NC_DOUBLE 1055 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini) 1056 #else 1057 ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini) 1058 #endif 1059 if (ierr.ne.0) then 1060 print *, NF_STRERROR(ierr) 1061 STOP 1062 ENDIF 1063 ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep) 1064 if (ierr.ne.0) then 1065 print *, NF_STRERROR(ierr) 1066 STOP 1067 ENDIF 1068 print*,'variable ', namedim, 'dimension ', lmdep 1069 ierr = NF_INQ_VARID(ncid,namedim,dimid) 1070 if (ierr.ne.0) then 1071 print *, NF_STRERROR(ierr) 1072 STOP 1073 ENDIF 1074 #ifdef NC_DOUBLE 1075 ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord) 1076 #else 1077 ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord) 1078 #endif 1079 if (ierr.ne.0) then 1080 print *, NF_STRERROR(ierr) 1081 STOP 1082 ENDIF 1083 c 1084 ALLOCATE ( champ(imdep*jmdep) ) 1085 1086 DO l = 1, lmdep 1087 dimfirst(1) = 1 1088 dimfirst(2) = 1 1089 dimfirst(3) = l 1090 c 1091 dimlast(1) = imdep 1092 dimlast(2) = jmdep 1093 dimlast(3) = 1 1094 c 1095 PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l) 1096 #ifdef NC_DOUBLE 1097 ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ) 1098 #else 1099 ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ) 1100 #endif 1101 if (ierr.ne.0) then 1102 print *, NF_STRERROR(ierr) 1103 STOP 1104 ENDIF 1105 1106 title='Albedo Amip' 1107 c 1108 CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini, 1109 . dlon, dlat, champ, interbar ) 1110 c 1111 c 1112 IF ( interbar ) THEN 1113 IF( l.EQ.1 ) THEN 1114 WRITE(6,*) '-------------------------------------------------', 1115 ,'------------------------' 1116 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 1117 , ' pour l Albedo Amip $$$ ' 1118 WRITE(6,*) '-------------------------------------------------', 1119 ,'------------------------' 1120 ENDIF 1121 1122 CALL inter_barxy ( imdep,jmdep -1,dlon, dlat , 1123 , champ, iim, jjm, rlonu, rlatv, jjp1, champint ) 1124 ELSE 1125 CALL grille_m(imdep, jmdep, dlon, dlat, champ, 1126 . iim, jjp1, rlonv, rlatu, champint ) 1127 ENDIF 1128 c 1129 DO j = 1,jjp1 1130 DO i = 1, iim 1131 champtime (i, j, l) = champint(i, j) 1132 ENDDO 1133 ENDDO 1134 ENDDO 1135 c 1136 DO l = 1, lmdep 1137 timeyear(l) = timecoord(l) 1138 ENDDO 1139 print 222, timeyear(:lmdep) 1140 c 1141 C interpolation temporelle 1142 DO j = 1, jjp1 1143 DO i = 1, iim 1144 DO l = 1, lmdep 1145 ax(l) = timeyear(l) 1146 ay(l) = champtime (i, j, l) 1147 ENDDO 1148 CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder) 1149 DO k = 1, 360 1150 time = FLOAT(k-1) 1151 CALL SPLINT(ax,ay,yder,lmdep,time,by) 1152 champan(i,j,k) = by 1153 ENDDO 1154 ENDDO 1155 ENDDO 1156 DO k = 1, 360 1157 DO j = 1, jjp1 1158 champan(iip1, j, k) = champan(1, j, k) 1159 ENDDO 1160 IF ( k.EQ.10 ) THEN 1161 DO j = 1, jjp1 1162 CALL minmax( iip1,champan(1,j,10),chmin,chmax ) 1163 PRINT *,' Albedo au temps 10 ', chmin,chmax,j 1164 ENDDO 1165 ENDIF 1166 ENDDO 1167 c 1168 DO k = 1, 360 1169 CALL gr_dyn_fi(1, iip1, jjp1, klon, 1170 . champan(1,1,k), phy_alb(1,k)) 1171 ENDDO 1172 c 1173 ierr = NF_CLOSE(ncid) 1174 c 1175 c 1176 DO k = 1, 360 1177 DO i = 1, klon 1178 phy_bil(i,k) = 0.0 1179 ENDDO 1180 ENDDO 1181 c 1182 PRINT*, 'Ecriture du fichier limit' 1183 c 1184 ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid) 1185 c 1186 ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30, 1187 . "Fichier conditions aux limites") 1188 ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim) 1189 ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim) 1190 c 1191 dims(1) = ndim 1192 dims(2) = ntim 1193 c 1194 #ifdef NC_DOUBLE 1195 ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim) 1196 #else 1197 ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim) 1198 #endif 1199 ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17, 1200 . "Jour dans l annee") 1201 IF (newlmt) THEN 1202 c 1203 #ifdef NC_DOUBLE 1204 ierr = NF_DEF_VAR (nid, "FOCE", NF_DOUBLE, 2,dims, id_FOCE) 1205 #else 1206 ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE) 1207 #endif 1208 ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14, 1209 . "Fraction ocean") 1210 c 1211 #ifdef NC_DOUBLE 1212 ierr = NF_DEF_VAR (nid, "FSIC", NF_DOUBLE, 2,dims, id_FSIC) 1213 #else 1214 ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC) 1215 #endif 1216 ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21, 1217 . "Fraction glace de mer") 1218 c 1219 #ifdef NC_DOUBLE 1220 ierr = NF_DEF_VAR (nid, "FTER", NF_DOUBLE, 2,dims, id_FTER) 1221 #else 1222 ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER) 1223 #endif 1224 ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14, 1225 . "Fraction terre") 1226 c 1227 #ifdef NC_DOUBLE 1228 ierr = NF_DEF_VAR (nid, "FLIC", NF_DOUBLE, 2,dims, id_FLIC) 1229 #else 1230 ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC) 1231 #endif 1232 ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17, 1233 . "Fraction land ice") 1234 c 1235 ELSE 1236 #ifdef NC_DOUBLE 1237 ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT) 1238 #else 1239 ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT) 1240 #endif 1241 ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23, 1242 . "Nature du sol (0,1,2,3)") 1243 ENDIF 1244 #ifdef NC_DOUBLE 1245 ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST) 1246 #else 1247 ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST) 1248 #endif 1249 ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35, 1250 . "Temperature superficielle de la mer") 1251 #ifdef NC_DOUBLE 1252 ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS) 1253 #else 1254 ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS) 1255 #endif 1256 ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32, 1257 . "Reference flux de chaleur au sol") 1258 #ifdef NC_DOUBLE 1259 ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB) 1260 #else 1261 ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB) 1262 #endif 1263 ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19, 1264 . "Albedo a la surface") 1265 #ifdef NC_DOUBLE 1266 ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG) 1267 #else 1268 ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG) 1269 #endif 1270 ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8, 1271 . "Rugosite") 1272 c 1273 ierr = NF_ENDDEF(nid) 1274 c 1275 DO k = 1, 360 1276 c 1277 debut(1) = 1 1278 debut(2) = k 1279 epais(1) = klon 1280 epais(2) = 1 1281 c 1282 #ifdef NC_DOUBLE 1283 ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k)) 1284 c 1285 IF (newlmt ) THEN 1286 ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais 1287 $ ,pctsrf_t(1,is_oce,k)) 1288 ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais 1289 $ ,pctsrf_t(1,is_sic,k)) 1290 ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais 1291 $ ,pctsrf_t(1,is_ter,k)) 1292 ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais 1293 $ ,pctsrf_t(1,is_lic,k)) 1294 ELSE 1295 ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais 1296 $ ,phy_nat(1,k)) 1297 ENDIF 1298 c 1299 ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k)) 1300 ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k)) 1301 ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k)) 1302 ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k)) 1303 #else 1304 ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k)) 1305 IF (newlmt ) THEN 1306 ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais 1307 $ ,pctsrf_t(1,is_oce,k)) 1308 ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais 1309 $ ,pctsrf_t(1,is_sic,k)) 1310 ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais 1311 $ ,pctsrf_t(1,is_ter,k)) 1312 ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais 1313 $ ,pctsrf_t(1,is_lic,k)) 1314 ELSE 1315 ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais 1316 $ ,phy_nat(1,k)) 1317 ENDIF 1318 ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k)) 1319 ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k)) 1320 ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k)) 1321 ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k)) 1322 #endif 1323 c 1324 ENDDO 1325 c 1326 ierr = NF_CLOSE(nid) 1327 c 1328 #else 1329 WRITE(lunout,*) 1330 & 'limit_netcdf: Earth-specific routine, needs Earth physics' 1331 #endif 1332 ! of #ifdef CPP_EARTH 1333 STOP 1334 END 306 !------------------------------------------------------------------------------- 307 ! Arguments: 308 CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name 309 CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB 310 LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels 311 INTEGER, INTENT(IN) :: ndays ! current year number of days 312 REAL, POINTER, DIMENSION(:,:) :: champo ! output field = f(t) 313 LOGICAL, OPTIONAL, INTENT(IN) :: flag 314 ! flag=T means: extrapolation (SST case) or old ice (SIC case) 315 REAL, OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask 316 !------------------------------------------------------------------------------- 317 ! Local variables: 318 !--- NetCDF 319 INTEGER :: ierr, ncid, varid ! NetCDF identifiers 320 CHARACTER(LEN=30) :: dnam ! dimension name 321 CHARACTER(LEN=80) :: varname ! NetCDF variable name 322 !--- dimensions 323 INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers 324 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini ! initial longitudes vector 325 REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini ! initial latitudes vector 326 REAL, POINTER, DIMENSION(:) :: dlon, dlat ! reordered lon/lat vectors 327 !--- fields 328 INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ' 329 REAL, ALLOCATABLE, DIMENSION(:) :: champ ! wanted field on initial grid 330 REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear 331 REAL, DIMENSION(iim,jjp1) :: champint ! interpolated field 332 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: champtime 333 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: champan 334 REAL :: time, by 335 !--- input files 336 CHARACTER(LEN=20) :: cal_in ! calendar 337 INTEGER :: ndays_in ! number of days 338 !--- misc 339 INTEGER :: i, j, k, l ! loop counters 340 REAL, ALLOCATABLE, DIMENSION(:,:) :: work ! used for extrapolation 341 CHARACTER(LEN=25) :: title ! for messages 342 LOGICAL :: extrp ! flag for extrapolation 343 REAL :: chmin, chmax 344 !------------------------------------------------------------------------------- 345 !---Variables depending on keyword 'mode' -------------------------------------- 346 NULLIFY(champo) 347 SELECT CASE(mode) 348 CASE('RUG'); varname='RUGOS'; title='Rugosite' 349 CASE('SIC'); varname='sicbcs'; title='Sea-ice' 350 CASE('SST'); varname='tosbcs'; title='SST' 351 CASE('ALB'); varname='ALBEDO'; title='Albedo' 352 END SELECT 353 extrp=(flag.AND.mode=='SST') 354 355 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------ 356 ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid); CALL ncerr(ierr,fnam) 357 ierr=NF90_INQ_VARID(ncid,varname,varid); CALL ncerr(ierr,fnam) 358 ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids); CALL ncerr(ierr,fnam) 359 360 !--- Longitude 361 ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep) 362 CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep)) 363 ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam) 364 ierr=NF90_GET_VAR(ncid,varid,dlon_ini); CALL ncerr(ierr,fnam) 365 WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep 366 367 !--- Latitude 368 ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep) 369 CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep)) 370 ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam) 371 ierr=NF90_GET_VAR(ncid,varid,dlat_ini); CALL ncerr(ierr,fnam) 372 WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep 373 374 !--- Time (variable is not needed - it is rebuilt - but calendar is) 375 ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep) 376 CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep)) 377 ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam) 378 cal_in=' ' 379 ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in) 380 IF(ierr/=NF90_NOERR) THEN 381 SELECT CASE(mode) 382 CASE('RUG','ALB'); cal_in='360d' 383 CASE('SIC','SST'); cal_in='gregorian' 384 END SELECT 385 WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d& 386 &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.' 387 END IF 388 WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in 389 390 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ---------------------- 391 !--- Determining input file number of days, depending on calendar 392 ndays_in=year_len(anneeref,cal_in) 393 394 !--- Time vector reconstruction (time vector from file is not trusted) 395 !--- If input records are not monthly, time sampling has to be constant ! 396 timeyear=mid_months(anneeref,cal_in,lmdep) 397 IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode) & 398 //' ne comportent pas 12, mais ',lmdep,' renregistrements.' 399 400 !--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------ 401 ALLOCATE(champ(imdep*jmdep),champtime(iim,jjp1,lmdep)) 402 IF(extrp) ALLOCATE(work(imdep,jmdep)) 403 404 WRITE(lunout,*) 405 WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.' 406 ierr=NF90_INQ_VARID(ncid,varname,varid); CALL ncerr(ierr,fnam) 407 DO l=1,lmdep 408 ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/)) 409 CALL ncerr(ierr,fnam) 410 CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar) 411 412 IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) 413 414 IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN 415 IF(l==1) THEN 416 WRITE(lunout,*) & 417 '-------------------------------------------------------------------------' 418 WRITE(lunout,*) & 419 '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' 420 WRITE(lunout,*) & 421 '-------------------------------------------------------------------------' 422 END IF 423 IF(mode=='RUG') champ=LOG(champ) 424 CALL inter_barxy(imdep, jmdep-1, dlon, dlat, champ, iim, jjm, rlonu, & 425 rlatv, jjp1, champint) 426 IF(mode=='RUG') THEN 427 champint=EXP(champint) 428 WHERE(NINT(mask)/=1) champint=0.001 429 END IF 430 ELSE 431 SELECT CASE(mode) 432 CASE('RUG'); CALL rugosite(imdep, jmdep, dlon, dlat, champ, & 433 iim, jjp1, rlonv, rlatu, champint, mask) 434 CASE('SIC'); CALL sea_ice (imdep, jmdep, dlon, dlat, champ, & 435 iim, jjp1, rlonv, rlatu, champint) 436 CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, & 437 iim, jjp1, rlonv, rlatu, champint) 438 END SELECT 439 END IF 440 champtime(:,:,l)=champint 441 END DO 442 ierr=NF90_CLOSE(ncid); CALL ncerr(ierr,fnam) 443 444 DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ) 445 IF(extrp) DEALLOCATE(work) 446 447 !--- TIME INTERPOLATION -------------------------------------------------------- 448 WRITE(lunout,*) 449 WRITE(lunout,*)'INTERPOLATION TEMPORELLE.' 450 WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear 451 WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays 452 ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays)) 453 DO j=1,jjp1 454 DO i=1,iim 455 CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder) 456 DO k=1,ndays 457 time=FLOAT((k-1)*ndays_in)/FLOAT(ndays) 458 CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by) 459 champan(i,j,k) = by 460 END DO 461 END DO 462 END DO 463 champan(iip1,:,:)=champan(1,:,:) 464 DEALLOCATE(yder,champtime,timeyear) 465 466 !--- Checking the result 467 DO j=1,jjp1 468 CALL minmax(iip1,champan(1,j,10),chmin,chmax) 469 WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j 470 END DO 471 472 !--- SPECIAL FILTER FOR SST: SST>271.38 ---------------------------------------- 473 IF(mode=='SST') THEN 474 WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38' 475 WHERE(champan<271.38) champan=271.38 476 END IF 477 478 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 --------------------------------------- 479 IF(mode=='SIC') THEN 480 WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 481 champan(:,:,:)=champan(:,:,:)/100. 482 champan(iip1,:,:)=champan(1,:,:) 483 WHERE(champan>1.0) champan=1.0 484 WHERE(champan<0.0) champan=0.0 485 END IF 486 487 write(*,*)'coin1' 488 !--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------ 489 ALLOCATE(champo(klon,ndays)) 490 DO k=1,ndays 491 CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k)) 492 END DO 493 write(*,*)'coin2' 494 DEALLOCATE(champan) 495 write(*,*)'coin3' 496 END SUBROUTINE get_2Dfield 497 ! 498 !------------------------------------------------------------------------------- 499 500 501 502 !------------------------------------------------------------------------------- 503 ! 504 FUNCTION year_len(y,cal_in) 505 ! 506 !------------------------------------------------------------------------------- 507 USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len 508 IMPLICIT NONE 509 !------------------------------------------------------------------------------- 510 ! Arguments: 511 INTEGER :: year_len 512 INTEGER, INTENT(IN) :: y 513 CHARACTER(LEN=*), INTENT(IN) :: cal_in 514 !------------------------------------------------------------------------------- 515 ! Local variables: 516 CHARACTER(LEN=20) :: cal_out ! calendar (for outputs) 517 !------------------------------------------------------------------------------- 518 !--- Getting the input calendar to reset at the end of the function 519 CALL ioget_calendar(cal_out) 520 521 !--- Unlocking calendar and setting it to wanted one 522 CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in)) 523 524 !--- Getting the number of days in this year 525 year_len=ioget_year_len(y) 526 527 !--- Back to original calendar 528 CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out)) 529 530 END FUNCTION year_len 531 ! 532 !------------------------------------------------------------------------------- 533 534 535 !------------------------------------------------------------------------------- 536 ! 537 FUNCTION mid_months(y,cal_in,nm) 538 ! 539 !------------------------------------------------------------------------------- 540 USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len 541 IMPLICIT NONE 542 !------------------------------------------------------------------------------- 543 ! Arguments: 544 INTEGER, INTENT(IN) :: y ! year 545 CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar 546 INTEGER, INTENT(IN) :: nm ! months/year number 547 REAL, DIMENSION(nm) :: mid_months ! mid-month times 548 !------------------------------------------------------------------------------- 549 ! Local variables: 550 CHARACTER(LEN=99) :: mess ! error message 551 CHARACTER(LEN=20) :: cal_out ! calendar (for outputs) 552 INTEGER, DIMENSION(nm) :: mnth ! months lengths (days) 553 INTEGER :: m ! months counter 554 INTEGER :: nd ! number of days 555 !------------------------------------------------------------------------------- 556 nd=year_len(y,cal_in) 557 558 IF(nm==12) THEN 559 560 !--- Getting the input calendar to reset at the end of the function 561 CALL ioget_calendar(cal_out) 562 563 !--- Unlocking calendar and setting it to wanted one 564 CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in)) 565 566 !--- Getting the length of each month 567 DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO 568 569 !--- Back to original calendar 570 CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out)) 571 572 ELSE IF(MODULO(nd,nm)/=0) THEN 573 WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',& 574 nm,' months/year. Months number should divide days number.' 575 CALL abort_gcm('mid_months',TRIM(mess),1) 576 577 ELSE 578 mnth=(/(m,m=1,nm,nd/nm)/) 579 END IF 580 581 !--- Mid-months times 582 mid_months(1)=0.5*FLOAT(mnth(1)) 583 DO k=2,nm 584 mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k)) 585 END DO 586 587 END FUNCTION mid_months 588 ! 589 !------------------------------------------------------------------------------- 590 591 592 !------------------------------------------------------------------------------- 593 ! 594 SUBROUTINE ncerr(ncres,fnam) 595 ! 596 !------------------------------------------------------------------------------- 597 ! Purpose: NetCDF errors handling. 598 !------------------------------------------------------------------------------- 599 USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR 600 IMPLICIT NONE 601 !------------------------------------------------------------------------------- 602 ! Arguments: 603 INTEGER, INTENT(IN) :: ncres 604 CHARACTER(LEN=*), INTENT(IN) :: fnam 605 !------------------------------------------------------------------------------- 606 #include "iniprint.h" 607 IF(ncres/=NF90_NOERR) THEN 608 WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.' 609 CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1) 610 END IF 611 612 END SUBROUTINE ncerr 613 ! 614 !------------------------------------------------------------------------------- 615 616 617 END SUBROUTINE limit_netcdf -
LMDZ4/trunk/libf/dyn3d/logic.h
r1046 r1319 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 4 ! … … 9 9 COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys, & 10 10 & statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus & 11 & ,read_start,ok_guide,ok_strato,ok_gradsfile 11 & ,read_start,ok_guide,ok_strato,ok_gradsfile & 12 & ,ok_limit,ok_etat0 12 13 13 14 LOGICAL purmats,forward,leapf,apphys,statcl,conser, & 14 15 & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus & 15 & ,read_start,ok_guide,ok_strato,ok_gradsfile 16 & ,read_start,ok_guide,ok_strato,ok_gradsfile & 17 & ,ok_limit,ok_etat0 16 18 17 19 INTEGER iflag_phys -
LMDZ4/trunk/libf/dyn3d/startvar.F90
r1318 r1319 2 2 ! $Id$ 3 3 ! 4 MODULE startvar 4 !******************************************************************************* 5 ! 6 MODULE startvar 7 ! 8 !******************************************************************************* 9 ! 10 !------------------------------------------------------------------------------- 11 ! Purpose: Access data from the database of atmospheric to initialize the model. 12 !------------------------------------------------------------------------------- 13 ! Comments: 14 ! 15 ! * This module is designed to work for Earth (and with ioipsl) 16 ! 17 ! * There are three ways to acces data, depending on the type of field 18 ! which needs to be extracted. In any case the call should come after a restget 19 ! and should be of the type : CALL startget(...) 20 ! 21 ! - A 1D variable on the physical grid : 22 ! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, & 23 ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) 24 ! 25 ! - A 2D variable on the dynamical grid : 26 ! CALL startget(varname, iml, jml, lon_in, lat_in, & 27 ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) 28 ! 29 ! - A 3D variable on the dynamical grid : 30 ! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, & 31 ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) 32 ! 33 ! * Data needs to be in NetCDF format 34 ! 35 ! * Variables should have the following names in the files: 36 ! 'RELIEF' : High resolution orography 37 ! 'ST' : Surface temperature 38 ! 'CDSW' : Soil moisture 39 ! 'Z' : Surface geopotential 40 ! 'SP' : Surface pressure 41 ! 'U' : East ward wind 42 ! 'V' : Northward wind 43 ! 'TEMP' : Temperature 44 ! 'R' : Relative humidity 45 ! 46 ! * There is a big mess with the longitude size. Should it be iml or iml+1 ? 47 ! I have chosen to use the iml+1 as an argument to this routine and we declare 48 ! internaly smaller fields when needed. This needs to be cleared once and for 49 ! all in LMDZ. A convention is required. 50 !------------------------------------------------------------------------------- 5 51 #ifdef CPP_EARTH 6 ! This module is designed to work for Earth (and with ioipsl) 7 ! 8 ! 9 ! There are three ways to access data from the database of atmospheric data which 10 ! can be used to initialize the model. This depends on the type of field which needs 11 ! to be extracted. In any case the call should come after a restget and should be of the type : 12 ! CALL startget(...) 13 ! 14 ! We will details the possible arguments to startget here : 15 ! 16 ! - A 2D variable on the dynamical grid : 17 ! CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar ) 18 ! 19 ! - A 1D variable on the physical grid : 20 ! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar ) 21 ! 22 ! 23 ! - A 3D variable on the dynamical grid : 24 ! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar ) 25 ! 26 ! 27 ! There is special constraint on the atmospheric data base except that the 28 ! the data needs to be in netCDF and the variables should have the the following 29 ! names in the file : 30 ! 31 ! 'RELIEF' : High resolution orography 32 ! 'ST' : Surface temperature 33 ! 'CDSW' : Soil moisture 34 ! 'Z' : Surface geopotential 35 ! 'SP' : Surface pressure 36 ! 'U' : East ward wind 37 ! 'V' : Northward wind 38 ! 'TEMP' : Temperature 39 ! 'R' : Relative humidity 40 ! 41 USE ioipsl 42 ! 43 ! 44 IMPLICIT NONE 45 ! 46 ! 47 PRIVATE 48 PUBLIC startget 49 ! 50 ! 51 INTERFACE startget 52 MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn 53 END INTERFACE 54 ! 55 INTEGER, SAVE :: fid_phys, fid_dyn 56 INTEGER, SAVE :: iml_phys, iml_rel, iml_dyn 57 INTEGER, SAVE :: jml_phys, jml_rel, jml_dyn 58 INTEGER, SAVE :: llm_dyn, ttm_dyn 59 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lon_phys, lon_rug, 60 . lon_alb, lon_rel, lon_dyn 61 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lat_phys, lat_rug, 62 . lat_alb, lat_rel, lat_dyn 63 REAL, ALLOCATABLE, SAVE, DIMENSION (:) :: levdyn_ini 64 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: relief, zstd, zsig, 65 . zgam, zthe, zpic, zval 66 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rugo, masque, phis 67 ! 68 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: tsol, qsol, psol_dyn 69 REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: var_ana3d 70 ! 71 CONTAINS 72 ! 73 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 74 ! 75 SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, 76 . champ, val_exp, jml2, lon_in2, lat_in2 , interbar, masque_lu ) 77 ! 78 ! There is a big mess with the size in logitude, should it be iml or iml+1. 79 ! I have chosen to use the iml+1 as an argument to this routine and we declare 80 ! internaly smaler fields when needed. This needs to be cleared once and for all in LMDZ. 81 ! A convention is required. 82 ! 83 ! 84 CHARACTER*(*), INTENT(in) :: varname 85 INTEGER, INTENT(in) :: iml, jml ,jml2 86 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 87 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 88 REAL, INTENT(inout) :: champ(iml,jml) 89 REAL, INTENT(in) :: val_exp 90 REAL, INTENT(in), optional :: masque_lu(iml,jml) 91 LOGICAL interbar 92 ! 93 ! This routine only works if the variable does not exist or is constant 94 ! 95 IF ( MINVAL(champ(:,:)).EQ.MAXVAL(champ(:,:)) .AND. 96 .MINVAL(champ(:,:)).EQ.val_exp ) THEN 97 ! 98 SELECTCASE(varname) 99 ! 100 CASE ('relief') 101 ! 102 ! If we do not have the orography we need to get it 103 ! 104 IF ( .NOT.ALLOCATED(relief)) THEN 105 ! 106 if (present(masque_lu)) then 107 CALL start_init_orog( iml, jml, lon_in, lat_in, 108 . jml2,lon_in2,lat_in2, interbar, masque_lu ) 109 else 110 CALL start_init_orog( iml, jml, lon_in, lat_in, 111 . jml2,lon_in2,lat_in2, interbar) 112 endif 113 ! 114 ENDIF 115 ! 116 IF (SIZE(relief) .NE. SIZE(champ)) THEN 117 ! 118 WRITE(*,*) 'STARTVAR module has been', 119 .' initialized to the wrong size' 120 STOP 121 ! 122 ENDIF 123 ! 124 champ(:,:) = relief(:,:) 125 ! 126 CASE ('rugosite') 127 ! 128 ! If we do not have the orography we need to get it 129 ! 130 IF ( .NOT.ALLOCATED(rugo)) THEN 131 ! 132 CALL start_init_orog( iml, jml, lon_in, lat_in, 133 . jml2,lon_in2,lat_in2 , interbar ) 134 ! 135 ENDIF 136 ! 137 IF (SIZE(rugo) .NE. SIZE(champ)) THEN 138 ! 139 WRITE(*,*) 140 . 'STARTVAR module has been initialized to the wrong size' 141 STOP 142 ! 143 ENDIF 144 ! 145 champ(:,:) = rugo(:,:) 146 ! 147 CASE ('masque') 148 ! 149 ! If we do not have the orography we need to get it 150 ! 151 IF ( .NOT.ALLOCATED(masque)) THEN 152 ! 153 CALL start_init_orog( iml, jml, lon_in, lat_in, 154 . jml2,lon_in2,lat_in2 , interbar ) 155 ! 156 ENDIF 157 ! 158 IF (SIZE(masque) .NE. SIZE(champ)) THEN 159 ! 160 WRITE(*,*) 161 . 'STARTVAR module has been initialized to the wrong size' 162 STOP 163 ! 164 ENDIF 165 ! 166 champ(:,:) = masque(:,:) 167 ! 168 CASE ('surfgeo') 169 ! 170 ! If we do not have the orography we need to get it 171 ! 172 IF ( .NOT.ALLOCATED(phis)) THEN 173 ! 174 CALL start_init_orog( iml, jml, lon_in, lat_in, 175 . jml2,lon_in2, lat_in2 , interbar ) 176 ! 177 ENDIF 178 ! 179 IF (SIZE(phis) .NE. SIZE(champ)) THEN 180 ! 181 WRITE(*,*) 182 . 'STARTVAR module has been initialized to the wrong size' 183 STOP 184 ! 185 ENDIF 186 ! 187 champ(:,:) = phis(:,:) 188 ! 189 CASE ('psol') 190 ! 191 ! If we do not have the orography we need to get it 192 ! 193 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 194 ! 195 CALL start_init_dyn( iml, jml, lon_in, lat_in, 196 . jml2,lon_in2, lat_in2 , interbar ) 197 ! 198 ENDIF 199 ! 200 IF (SIZE(psol_dyn) .NE. SIZE(champ)) THEN 201 ! 202 WRITE(*,*) 203 . 'STARTVAR module has been initialized to the wrong size' 204 STOP 205 ! 206 ENDIF 207 ! 208 champ(:,:) = psol_dyn(:,:) 209 ! 210 CASE DEFAULT 211 ! 212 WRITE(*,*) 'startget_phys2d' 213 WRITE(*,*) 'No rule is present to extract variable', 214 . varname(:LEN_TRIM(varname)),' from any data set' 215 STOP 216 ! 217 END SELECT 218 ! 219 ELSE 220 ! 221 ! There are a few fields we might need if we need to interpolate 3D filed. Thus if they come through here we 222 ! will catch them 223 ! 224 SELECTCASE(varname) 225 ! 226 CASE ('surfgeo') 227 ! 228 IF ( .NOT.ALLOCATED(phis)) THEN 229 ALLOCATE(phis(iml,jml)) 230 ENDIF 231 ! 232 IF (SIZE(phis) .NE. SIZE(champ)) THEN 233 ! 234 WRITE(*,*) 235 . 'STARTVAR module has been initialized to the wrong size' 236 STOP 237 ! 238 ENDIF 239 ! 240 phis(:,:) = champ(:,:) 241 ! 242 END SELECT 243 ! 244 ENDIF 245 ! 246 END SUBROUTINE startget_phys2d 247 ! 248 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 249 ! 250 SUBROUTINE start_init_orog ( iml,jml,lon_in, lat_in,jml2,lon_in2 , 251 , lat_in2 , interbar, masque_lu ) 252 ! 253 INTEGER, INTENT(in) :: iml, jml, jml2 254 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 255 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 256 REAL, intent(in), optional :: masque_lu(iml,jml) 257 LOGICAL interbar 258 ! 259 ! LOCAL 260 ! 261 LOGICAL interbar2 262 REAL :: lev(1), date, dt,chmin,chmax 263 INTEGER :: itau(1), fid 264 INTEGER :: llm_tmp, ttm_tmp 265 INTEGER :: i, j 266 INTEGER :: iret 267 CHARACTER*25 title 268 REAL, ALLOCATABLE :: relief_hi(:,:) 269 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 270 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) 271 REAL, ALLOCATABLE :: tmp_var(:,:) 272 INTEGER, ALLOCATABLE :: tmp_int(:,:) 273 ! 274 CHARACTER*120 :: orogfname 275 LOGICAL :: check=.TRUE. 276 ! 277 ! 278 orogfname = 'Relief.nc' 279 ! 280 IF ( check ) WRITE(*,*) 'Reading the high resolution orography' 281 ! 282 CALL flininfo(orogfname,iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 283 ! 284 ALLOCATE (lat_rel(iml_rel,jml_rel), stat=iret) 285 ALLOCATE (lon_rel(iml_rel,jml_rel), stat=iret) 286 ALLOCATE (relief_hi(iml_rel,jml_rel), stat=iret) 287 ! 288 CALL flinopen(orogfname, .FALSE., iml_rel, jml_rel, 289 .llm_tmp, lon_rel, lat_rel, lev, ttm_tmp, 290 . itau, date, dt, fid) 291 ! 292 CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp, 293 . ttm_tmp, 1, 1, relief_hi) 294 ! 295 CALL flinclo(fid) 296 ! 297 ! In case we have a file which is in degrees we do the transformation 298 ! 299 ALLOCATE(lon_rad(iml_rel)) 300 ALLOCATE(lon_ini(iml_rel)) 301 302 IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 303 lon_ini(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0 304 ELSE 305 lon_ini(:) = lon_rel(:,1) 306 ENDIF 307 308 ALLOCATE(lat_rad(jml_rel)) 309 ALLOCATE(lat_ini(jml_rel)) 310 311 IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 312 lat_ini(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0 313 ELSE 314 lat_ini(:) = lat_rel(1,:) 315 ENDIF 316 ! 317 ! 318 319 title='RELIEF' 320 321 interbar2 = .FALSE. 322 CALL conf_dat2d(title,iml_rel, jml_rel, lon_ini, lat_ini, 323 . lon_rad, lat_rad, relief_hi , interbar2 ) 324 325 IF ( check ) WRITE(*,*) 'Computes all the parameters needed', 326 .' for the gravity wave drag code' 327 ! 328 ! Allocate the data we need to put in the interpolated fields 329 ! 330 ! RELIEF: orographie moyenne 331 ALLOCATE(relief(iml,jml)) 332 ! zphi : orographie moyenne 333 ALLOCATE(phis(iml,jml)) 334 ! zstd: deviation standard de l'orographie sous-maille 335 ALLOCATE(zstd(iml,jml)) 336 ! zsig: pente de l'orographie sous-maille 337 ALLOCATE(zsig(iml,jml)) 338 ! zgam: anisotropy de l'orographie sous maille 339 ALLOCATE(zgam(iml,jml)) 340 ! zthe: orientation de l'axe oriente dans la direction 341 ! de plus grande pente de l'orographie sous maille 342 ALLOCATE(zthe(iml,jml)) 343 ! zpic: hauteur pics de la SSO 344 ALLOCATE(zpic(iml,jml)) 345 ! zval: hauteur vallees de la SSO 346 ALLOCATE(zval(iml,jml)) 347 ! masque : Masque terre ocean 348 ALLOCATE(tmp_int(iml,jml)) 349 ALLOCATE(masque(iml,jml)) 350 351 masque = -99999. 352 if (present(masque_lu)) then 353 masque = masque_lu 354 endif 355 ! 356 CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, 357 . iml-1, jml, lon_in, lat_in, 358 . phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque) 359 phis = phis * 9.81 360 ! 361 ! masque(:,:) = FLOAT(tmp_int(:,:)) 362 ! 363 ! Compute surface roughness 364 ! 365 IF ( check ) WRITE(*,*) 366 .'Compute surface roughness induced by the orography' 367 ! 368 ALLOCATE(rugo(iml,jml)) 369 ALLOCATE(tmp_var(iml-1,jml)) 370 ! 371 CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, 372 . iml-1, jml, lon_in, lat_in, tmp_var) 373 ! 374 DO j = 1, jml 375 DO i = 1, iml-1 376 rugo(i,j) = tmp_var(i,j) 377 ENDDO 378 rugo(iml,j) = tmp_var(1,j) 379 ENDDO 380 c 381 cc *** rugo n'est pas utilise pour l'instant ****** 382 ! 383 ! Build land-sea mask 384 ! 385 ! 386 RETURN 387 ! 388 END SUBROUTINE start_init_orog 389 ! 390 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 391 ! 392 SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, 393 .lat_in, nbindex, champ, val_exp ,jml2, lon_in2, lat_in2,interbar) 394 ! 395 CHARACTER*(*), INTENT(in) :: varname 396 INTEGER, INTENT(in) :: iml, jml, nbindex, jml2 397 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 398 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 399 REAL, INTENT(inout) :: champ(nbindex) 400 REAL, INTENT(in) :: val_exp 401 LOGICAL interbar 402 ! 403 ! 404 ! This routine only works if the variable does not exist or is constant 405 ! 406 IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND. 407 .MINVAL(champ(:)).EQ.val_exp ) THEN 408 SELECTCASE(varname) 409 CASE ('tsol') 410 IF ( .NOT.ALLOCATED(tsol)) THEN 411 CALL start_init_phys( iml, jml, lon_in, lat_in, 412 . jml2, lon_in2, lat_in2, interbar ) 413 ENDIF 414 IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 415 WRITE(*,*) 416 . 'STARTVAR module has been initialized to the wrong size' 417 STOP 418 ENDIF 419 CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ) 420 CASE ('qsol') 421 IF ( .NOT.ALLOCATED(qsol)) THEN 422 CALL start_init_phys( iml, jml, lon_in, lat_in, 423 . jml2, lon_in2,lat_in2 , interbar ) 424 ENDIF 425 IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 426 WRITE(*,*) 427 . 'STARTVAR module has been initialized to the wrong size' 428 STOP 429 ENDIF 430 CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ) 431 CASE ('psol') 432 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 433 CALL start_init_dyn( iml, jml, lon_in, lat_in, 434 . jml2, lon_in2,lat_in2 , interbar ) 435 ENDIF 436 IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN 437 WRITE(*,*) 438 . 'STARTVAR module has been initialized to the wrong size' 439 STOP 440 ENDIF 441 CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ) 442 CASE ('zmea') 443 IF ( .NOT.ALLOCATED(relief)) THEN 444 CALL start_init_orog( iml, jml, lon_in, lat_in, 445 . jml2, lon_in2,lat_in2 , interbar ) 446 ENDIF 447 IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 448 WRITE(*,*) 449 . 'STARTVAR module has been initialized to the wrong size' 450 STOP 451 ENDIF 452 CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ) 453 CASE ('zstd') 454 IF ( .NOT.ALLOCATED(zstd)) THEN 455 CALL start_init_orog( iml, jml, lon_in, lat_in, 456 . jml2, lon_in2,lat_in2 , interbar ) 457 ENDIF 458 IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 459 WRITE(*,*) 460 . 'STARTVAR module has been initialized to the wrong size' 461 STOP 462 ENDIF 463 CALL gr_dyn_fi(1, iml, jml, nbindex,zstd, champ) 464 CASE ('zsig') 465 IF ( .NOT.ALLOCATED(zsig)) THEN 466 CALL start_init_orog( iml, jml, lon_in, lat_in, 467 . jml2, lon_in2,lat_in2 , interbar ) 468 ENDIF 469 IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 470 WRITE(*,*) 471 . 'STARTVAR module has been initialized to the wrong size' 472 STOP 473 ENDIF 474 CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ) 475 CASE ('zgam') 476 IF ( .NOT.ALLOCATED(zgam)) THEN 477 CALL start_init_orog( iml, jml, lon_in, lat_in, 478 . jml2, lon_in2,lat_in2 , interbar ) 479 ENDIF 480 IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 481 WRITE(*,*) 482 . 'STARTVAR module has been initialized to the wrong size' 483 STOP 484 ENDIF 485 CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ) 486 CASE ('zthe') 487 IF ( .NOT.ALLOCATED(zthe)) THEN 488 CALL start_init_orog( iml, jml, lon_in, lat_in, 489 . jml2, lon_in2,lat_in2 , interbar ) 490 ENDIF 491 IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 492 WRITE(*,*) 493 . 'STARTVAR module has been initialized to the wrong size' 494 STOP 495 ENDIF 496 CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ) 497 CASE ('zpic') 498 IF ( .NOT.ALLOCATED(zpic)) THEN 499 CALL start_init_orog( iml, jml, lon_in, lat_in, 500 . jml2, lon_in2,lat_in2 , interbar ) 501 ENDIF 502 IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 503 WRITE(*,*) 504 . 'STARTVAR module has been initialized to the wrong size' 505 STOP 506 ENDIF 507 CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ) 508 CASE ('zval') 509 IF ( .NOT.ALLOCATED(zval)) THEN 510 CALL start_init_orog( iml, jml, lon_in, lat_in, 511 . jml2, lon_in2,lat_in2 , interbar ) 512 ENDIF 513 IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN 514 WRITE(*,*) 515 . 'STARTVAR module has been initialized to the wrong size' 516 STOP 517 ENDIF 518 CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ) 519 CASE ('rads') 520 champ(:) = 0.0 521 CASE ('snow') 522 champ(:) = 0.0 523 cIM "slab" ocean 524 CASE ('tslab') 525 champ(:) = 0.0 526 CASE ('seaice') 527 champ(:) = 0.0 528 CASE ('rugmer') 529 champ(:) = 0.001 530 CASE ('agsno') 531 champ(:) = 50.0 532 CASE DEFAULT 533 WRITE(*,*) 'startget_phys1d' 534 WRITE(*,*) 'No rule is present to extract variable ', 535 . varname(:LEN_TRIM(varname)),' from any data set' 536 STOP 537 END SELECT 538 ELSE 539 ! 540 ! If we see tsol we catch it as we may need it for a 3D interpolation 541 ! 542 SELECTCASE(varname) 543 CASE ('tsol') 544 IF ( .NOT.ALLOCATED(tsol)) THEN 545 ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) )) 546 ENDIF 547 CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol) 52 USE ioipsl 53 IMPLICIT NONE 54 55 PRIVATE 56 PUBLIC startget 57 INTERFACE startget 58 MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn 59 END INTERFACE 60 61 REAL, SAVE :: deg2rad, pi 62 INTEGER, SAVE :: iml_rel, jml_rel 63 INTEGER, SAVE :: fid_phys, iml_phys, jml_phys 64 INTEGER, SAVE :: fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn 65 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn 66 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn 67 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel 68 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel 69 REAL, DIMENSION(:), ALLOCATABLE, TARGET, SAVE :: levdyn_ini 70 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam 71 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval 72 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol 73 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: psol_dyn 74 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d 75 76 CONTAINS 77 78 !------------------------------------------------------------------------------- 79 ! 80 SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, & 81 val_exp ,jml2, lon_in2, lat_in2, ibar) 82 ! 83 !------------------------------------------------------------------------------- 84 ! Comment: 85 ! This routine only works if the variable does not exist or is constant. 86 !------------------------------------------------------------------------------- 87 ! Arguments: 88 CHARACTER(LEN=*), INTENT(IN) :: varname 89 INTEGER, INTENT(IN) :: iml, jml 90 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 91 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 92 INTEGER, INTENT(IN) :: nbindex 93 REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ 94 REAL, INTENT(IN) :: val_exp 95 INTEGER, INTENT(IN) :: jml2 96 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 97 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 98 LOGICAL, INTENT(IN) :: ibar 99 !------------------------------------------------------------------------------- 100 ! Local variables: 101 #include "iniprint.h" 102 REAL, DIMENSION(:,:), POINTER :: v2d 103 !------------------------------------------------------------------------------- 104 v2d=>NULL() 105 IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN 106 107 !--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES 108 SELECT CASE(varname) 109 CASE('tsol') 110 IF(.NOT.ALLOCATED(tsol)) & 111 CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 112 CASE('qsol') 113 IF(.NOT.ALLOCATED(qsol)) & 114 CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 115 CASE('psol') 116 IF(.NOT.ALLOCATED(psol_dyn)) & 117 CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 118 CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval') 119 IF(.NOT.ALLOCATED(relief)) & 120 CALL start_init_orog(iml,jml,lon_in,lat_in) 121 CASE('rads','snow','tslab','seaice','rugmer','agsno') 122 CASE DEFAULT 123 WRITE(lunout,*)'startget_phys1d' 124 WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& 125 //' from any data set'; STOP 126 END SELECT 127 128 !--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE 129 SELECT CASE(varname) 130 CASE('rads','snow','tslab','seaice'); champ=0.0 131 CASE('rugmer'); champ(:)=0.001 132 CASE('agsno'); champ(:)=50.0 133 CASE DEFAULT 134 SELECT CASE(varname) 135 CASE('tsol'); v2d=>tsol 136 CASE('qsol'); v2d=>qsol 137 CASE('psol'); v2d=>psol_dyn 138 CASE('zmea'); v2d=>relief 139 CASE('zstd'); v2d=>zstd 140 CASE('zsig'); v2d=>zsig 141 CASE('zgam'); v2d=>zgam 142 CASE('zthe'); v2d=>zthe 143 CASE('zpic'); v2d=>zpic 144 CASE('zval'); v2d=>zval 548 145 END SELECT 549 ENDIF 550 END SUBROUTINE startget_phys1d 551 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 552 ! 553 SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2, 554 . lon_in2, lat_in2 , interbar ) 555 ! 556 INTEGER, INTENT(in) :: iml, jml ,jml2 557 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 558 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 559 LOGICAL interbar 560 ! 561 ! LOCAL 562 ! 563 !ac REAL :: lev(1), date, dt 564 REAL :: date, dt 565 REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini 566 !ac 567 INTEGER :: itau(1) 568 INTEGER :: llm_tmp, ttm_tmp 569 INTEGER :: i, j 570 ! 571 CHARACTER*25 title 572 CHARACTER*120 :: physfname 573 LOGICAL :: check=.TRUE. 574 ! 575 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 576 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) 577 REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:) 578 ! 579 physfname = 'ECPHY.nc' 580 ! 581 IF ( check ) WRITE(*,*) 'Opening the surface analysis' 582 ! 583 CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, 584 . ttm_tmp, fid_phys) 585 ! 586 ALLOCATE (lat_phys(iml_phys,jml_phys)) 587 ALLOCATE (lon_phys(iml_phys,jml_phys)) 588 !ac 589 ALLOCATE (levphys_ini(llm_tmp)) 590 ! 591 ! CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 592 ! . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp, 593 ! . itau, date, dt, fid_phys) 594 ! 595 CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 596 . llm_tmp, lon_phys, lat_phys, levphys_ini, ttm_tmp, 597 . itau, date, dt, fid_phys) 598 ! 599 DEALLOCATE (levphys_ini) 600 !ac 601 ! 602 ! Allocate the space we will need to get the data out of this file 603 ! 604 ALLOCATE(var_ana(iml_phys, jml_phys)) 605 ! 606 ! In case we have a file which is in degrees we do the transformation 607 ! 608 ALLOCATE(lon_rad(iml_phys)) 609 ALLOCATE(lon_ini(iml_phys)) 610 611 IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 612 lon_ini(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0 613 ELSE 614 lon_ini(:) = lon_phys(:,1) 615 ENDIF 616 617 ALLOCATE(lat_rad(jml_phys)) 618 ALLOCATE(lat_ini(jml_phys)) 619 620 IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 621 lat_ini(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0 622 ELSE 623 lat_ini(:) = lat_phys(1,:) 624 ENDIF 625 626 627 ! 628 ! We get the two standard varibales 629 ! Surface temperature 630 ! 631 ALLOCATE(tsol(iml,jml)) 632 ALLOCATE(tmp_var(iml-1,jml)) 633 ! 634 ! 635 636 CALL flinget(fid_phys, 'ST', iml_phys, jml_phys, 637 .llm_tmp, ttm_tmp, 1, 1, var_ana) 638 639 title='ST' 640 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, 641 . lon_rad, lat_rad, var_ana , interbar ) 642 643 IF ( interbar ) THEN 644 WRITE(6,*) '-------------------------------------------------', 645 ,'--------------' 646 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 647 , ' pour ST $$$ ' 648 WRITE(6,*) '-------------------------------------------------', 649 ,'--------------' 650 CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad , 651 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var ) 652 ELSE 653 CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad, 654 . var_ana, iml-1, jml, lon_in, lat_in, tmp_var ) 655 ENDIF 656 657 CALL gr_int_dyn(tmp_var, tsol, iml-1, jml) 658 ! 659 ! Soil moisture 660 ! 661 ALLOCATE(qsol(iml,jml)) 662 CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys, 663 . llm_tmp, ttm_tmp, 1, 1, var_ana) 664 665 title='CDSW' 666 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, 667 . lon_rad, lat_rad, var_ana, interbar ) 668 669 IF ( interbar ) THEN 670 WRITE(6,*) '-------------------------------------------------', 671 ,'--------------' 672 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 673 , ' pour CDSW $$$ ' 674 WRITE(6,*) '-------------------------------------------------', 675 ,'--------------' 676 CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad , 677 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var ) 678 ELSE 679 CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad, 680 . var_ana, iml-1, jml, lon_in, lat_in, tmp_var ) 681 ENDIF 682 c 683 CALL gr_int_dyn(tmp_var, qsol, iml-1, jml) 684 ! 685 CALL flinclo(fid_phys) 686 ! 687 END SUBROUTINE start_init_phys 688 ! 689 ! 690 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 691 ! 692 ! 693 SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in, 694 . lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 , 695 , interbar ) 696 ! 697 ! ARGUMENTS 698 ! 699 CHARACTER*(*), INTENT(in) :: varname 700 INTEGER, INTENT(in) :: iml, jml, lml, jml2 701 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 702 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 703 REAL, INTENT(in) :: pls(iml, jml, lml) 704 REAL, INTENT(in) :: workvar(iml, jml, lml) 705 REAL, INTENT(inout) :: champ(iml, jml, lml) 706 REAL, INTENT(in) :: val_exp 707 LOGICAL interbar 708 ! 709 ! LOCAL 710 ! 711 INTEGER :: il, ij, ii 712 REAL :: xppn, xpps 713 ! 146 IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN 147 WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size' 148 STOP 149 END IF 150 CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ) 151 END SELECT 152 153 ELSE 154 155 !--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION 156 SELECT CASE(varname) 157 CASE('tsol') 158 IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml)) 159 CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol) 160 END SELECT 161 162 END IF 163 164 END SUBROUTINE startget_phys1d 165 ! 166 !------------------------------------------------------------------------------- 167 168 169 !------------------------------------------------------------------------------- 170 ! 171 SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp, & 172 jml2, lon_in2, lat_in2 , ibar, msk) 173 ! 174 !------------------------------------------------------------------------------- 175 ! Comment: 176 ! This routine only works if the variable does not exist or is constant. 177 !------------------------------------------------------------------------------- 178 ! Arguments: 179 CHARACTER(LEN=*), INTENT(IN) :: varname 180 INTEGER, INTENT(IN) :: iml, jml 181 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 182 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 183 REAL, DIMENSION(iml,jml), INTENT(INOUT) :: champ 184 REAL, INTENT(IN) :: val_exp 185 INTEGER, INTENT(IN) :: jml2 186 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 187 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 188 LOGICAL, INTENT(IN) :: ibar 189 REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk 190 !------------------------------------------------------------------------------- 191 ! Local variables: 192 #include "iniprint.h" 193 REAL, DIMENSION(:,:), POINTER :: v2d=>NULL() 194 LOGICAL :: lrelief1, lrelief2 195 !------------------------------------------------------------------------------- 196 v2d=>NULL() 197 lrelief1=(.NOT.ALLOCATED(relief).AND. PRESENT(msk)) 198 lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk)) 199 IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN 200 201 !--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES 202 SELECT CASE(varname) 203 CASE('psol') 204 IF(.NOT.ALLOCATED(psol_dyn)) & 205 CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 206 CASE('relief') 207 IF(lrelief1) CALL start_init_orog(iml,jml,lon_in,lat_in,msk) 208 IF(lrelief2) CALL start_init_orog(iml,jml,lon_in,lat_in) 209 CASE('surfgeo') 210 IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in) 211 CASE('rugosite','masque') 212 IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in) 213 CASE DEFAULT 214 WRITE(lunout,*)'startget_phys2d' 215 WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& 216 //' from any data set'; STOP 217 END SELECT 218 219 !--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE 220 SELECT CASE(varname) 221 CASE('psol'); v2d=>psol_dyn 222 CASE('relief'); v2d=>relief 223 CASE('rugosite'); v2d=>rugo 224 CASE('masque'); v2d=>masque 225 CASE('surfgeo'); v2d=>phis 226 END SELECT 227 IF(SIZE(champ)/=SIZE(v2d)) THEN 228 WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size' 229 STOP 230 END IF 231 232 champ(:,:)=v2d(:,:) 233 234 ELSE 235 236 !--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION 237 SELECT CASE(varname) 238 CASE ('surfgeo') 239 IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml)) 240 IF(SIZE(phis)/=SIZE(champ)) THEN 241 WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size' 242 STOP 243 END IF 244 phis(:,:)=champ(:,:) 245 END SELECT 246 247 END IF 248 249 END SUBROUTINE startget_phys2d 250 ! 251 !------------------------------------------------------------------------------- 252 253 254 !------------------------------------------------------------------------------- 255 ! 256 SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in, lml, pls,workvar,& 257 champ, val_exp, jml2, lon_in2, lat_in2, ibar) 258 !------------------------------------------------------------------------------- 259 ! Comment: 260 ! This routine only works if the variable does not exist or is constant. 261 !------------------------------------------------------------------------------- 262 ! Arguments: 263 CHARACTER(LEN=*), INTENT(IN) :: varname 264 INTEGER, INTENT(IN) :: iml, jml 265 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 266 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 267 INTEGER, INTENT(IN) :: lml 268 REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls, workvar 269 REAL, DIMENSION(iml,jml,lml), INTENT(INOUT) :: champ 270 REAL, INTENT(IN) :: val_exp 271 INTEGER, INTENT(IN) :: jml2 272 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 273 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 274 LOGICAL, INTENT(IN) :: ibar 275 !------------------------------------------------------------------------------- 276 ! Local variables: 277 #include "iniprint.h" 278 #include "dimensions.h" 279 #include "comconst.h" 280 #include "paramet.h" 281 #include "comgeom2.h" 282 REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL() 283 CHARACTER(LEN=10) :: vname 284 INTEGER :: il 285 REAL :: xppn, xpps 286 !------------------------------------------------------------------------------- 287 NULLIFY(v3d) 288 IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN 289 290 !--- READING UNALLOCATED FILES 291 IF(.NOT.ALLOCATED(psol_dyn)) & 292 CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 293 294 !--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS 295 SELECT CASE(varname) 296 CASE('u'); vname='U' 297 CASE('v'); vname='V' 298 CASE('t','tpot'); vname='TEMP' 299 CASE('q'); vname='R' 300 CASE DEFAULT 301 WRITE(lunout,*)'startget_dyn' 302 WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& 303 //' from any data set'; STOP 304 END SELECT 305 CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2, & 306 lon_in2, lat_in2, pls, champ,ibar ) 307 308 !--- COMPUTING THE REQUIRED FILED 309 SELECT CASE(varname) 310 CASE('u') !--- Eastward wind 311 DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO 312 champ(iml,:,:)=champ(1,:,:) 313 314 CASE('v') !--- Northward wind 315 DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO 316 champ(iml,:,:)=champ(1,:,:) 317 318 CASE('tpot') !--- Temperature 319 IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN 320 champ=champ*cpp/workvar 321 DO il=1,lml 322 xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln 323 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 324 champ(:,1 ,il) = xppn 325 champ(:,jml,il) = xpps 326 END DO 327 ELSE 328 WRITE(lunout,*)'Could not compute potential temperature as the' 329 WRITE(lunout,*)'Exner function is missing or constant.'; STOP 330 END IF 331 332 CASE('q') !--- Relat. humidity 333 IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN 334 champ=0.01*champ*workvar 335 WHERE(champ<0.) champ=1.0E-10 336 DO il=1,lml 337 xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln 338 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 339 champ(:,1 ,il) = xppn 340 champ(:,jml,il) = xpps 341 END DO 342 ELSE 343 WRITE(lunout,*)'Could not compute specific humidity as the' 344 WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP 345 END IF 346 347 END SELECT 348 349 END IF 350 351 END SUBROUTINE startget_dyn 352 ! 353 !------------------------------------------------------------------------------- 354 355 356 !------------------------------------------------------------------------------- 357 ! 358 SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu) 359 ! 360 !------------------------------------------------------------------------------- 361 ! Arguments: 362 INTEGER, INTENT(IN) :: iml, jml 363 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 364 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 365 REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu 366 !------------------------------------------------------------------------------- 367 ! Local variables: 368 #include "iniprint.h" 369 CHARACTER(LEN=25) :: title 370 CHARACTER(LEN=120) :: orofname 371 LOGICAL :: check=.TRUE. 372 REAL, DIMENSION(1) :: lev 373 REAL :: date, dt 374 INTEGER, DIMENSION(1) :: itau 375 INTEGER :: fid, llm_tmp, ttm_tmp 376 REAL, DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var 377 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 378 !------------------------------------------------------------------------------- 379 pi=2.0*ASIN(1.0); deg2rad=pi/180.0 380 381 orofname = 'Relief.nc'; title='RELIEF' 382 IF(check) WRITE(lunout,*)'Reading the high resolution orography' 383 CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 384 385 ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) 386 CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& 387 lev, ttm_tmp, itau, date, dt, fid) 388 ALLOCATE(relief_hi(iml_rel,jml_rel)) 389 CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) 390 CALL flinclo(fid) 391 392 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 393 ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) 394 lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad 395 lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad 396 397 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 398 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 399 CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, & 400 relief_hi, .FALSE.) 401 DEALLOCATE(lon_ini,lat_ini) 402 403 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro 404 IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra& 405 &g code' 406 407 ALLOCATE(phis(iml,jml)) ! Geopotentiel au sol 408 ALLOCATE(zstd(iml,jml)) ! Deviation standard de l'orographie sous-maille 409 ALLOCATE(zsig(iml,jml)) ! Pente de l'orographie sous-maille 410 ALLOCATE(zgam(iml,jml)) ! Anisotropie de l'orographie sous maille 411 ALLOCATE(zthe(iml,jml)) ! Orientation axe +grande pente d'oro sous maille 412 ALLOCATE(zpic(iml,jml)) ! Hauteur pics de la SSO 413 ALLOCATE(zval(iml,jml)) ! Hauteur vallees de la SSO 414 ALLOCATE(relief(iml,jml)) ! Orographie moyenne 415 ALLOCATE(masque(iml,jml)) ! Masque terre ocean 416 masque = -99999. 417 IF(PRESENT(masque_lu)) masque=masque_lu 418 419 CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, & 420 lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque) 421 phis = phis * 9.81 422 423 !--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! ) 424 IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography' 425 ALLOCATE(rugo (iml ,jml)) 426 ALLOCATE(tmp_var(iml-1,jml)) 427 CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, & 428 lon_in, lat_in, tmp_var) 429 rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:) 430 DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad) 431 RETURN 432 433 END SUBROUTINE start_init_orog 434 ! 435 !------------------------------------------------------------------------------- 436 437 438 !------------------------------------------------------------------------------- 439 ! 440 SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 441 ! 442 !------------------------------------------------------------------------------- 443 ! Arguments: 444 INTEGER, INTENT(IN) :: iml, jml 445 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 446 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 447 INTEGER, INTENT(IN) :: jml2 448 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 449 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 450 LOGICAL, INTENT(IN) :: ibar 451 !------------------------------------------------------------------------------- 452 ! Local variables: 453 #include "iniprint.h" 454 CHARACTER(LEN=25) :: title 455 CHARACTER(LEN=120) :: physfname 456 LOGICAL :: check=.TRUE. 457 REAL :: date, dt 458 INTEGER, DIMENSION(1) :: itau 459 INTEGER :: llm_tmp, ttm_tmp 460 REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana 461 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 462 REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini 463 !------------------------------------------------------------------------------- 464 physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0 465 IF(check) WRITE(lunout,*)'Opening the surface analysis' 466 CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys) 467 468 ALLOCATE(lat_phys(iml_phys,jml_phys)) 469 ALLOCATE(lon_phys(iml_phys,jml_phys)) 470 ALLOCATE(levphys_ini(llm_tmp)) 471 CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys, & 472 lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys) 473 DEALLOCATE(levphys_ini) 474 475 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 476 ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys)) 477 lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad 478 lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad 479 480 ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) 481 482 !--- SURFACE TEMPERATURE 483 title='ST' 484 ALLOCATE(tsol(iml,jml)) 485 CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana) 486 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,& 487 var_ana , ibar ) 488 CALL interp_startvar(title, ibar, .TRUE., & 489 iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 490 lon_in, lat_in, lon_in2, lat_in2, tsol) 491 492 !--- SOIL MOISTURE 493 title='CDSW' 494 ALLOCATE(qsol(iml,jml)) 495 CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana) 496 CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,& 497 var_ana, ibar ) 498 CALL interp_startvar(title, ibar, .TRUE., & 499 iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 500 lon_in, lat_in, lon_in2, lat_in2, qsol) 501 502 CALL flinclo(fid_phys) 503 504 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) 505 506 END SUBROUTINE start_init_phys 507 ! 508 !------------------------------------------------------------------------------- 509 510 511 !------------------------------------------------------------------------------- 512 ! 513 SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 514 ! 515 !------------------------------------------------------------------------------- 516 ! Arguments: 517 INTEGER, INTENT(IN) :: iml, jml 518 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 519 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 520 INTEGER, INTENT(IN) :: jml2 521 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 522 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 523 LOGICAL, INTENT(IN) :: ibar 524 !------------------------------------------------------------------------------- 525 ! Local variables: 526 #include "iniprint.h" 714 527 #include "dimensions.h" 715 528 #include "paramet.h" 716 529 #include "comgeom2.h" 717 #include "comconst.h" 718 ! 719 ! This routine only works if the variable does not exist or is constant 720 ! 721 IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND. 722 . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN 723 ! 724 SELECTCASE(varname) 725 CASE ('u') 726 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 727 CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 , 728 . lon_in2,lat_in2 , interbar ) 729 ENDIF 730 CALL start_inter_3d('U', iml, jml, lml, lon_in, 731 . lat_in, jml2, lon_in2, lat_in2, pls, champ,interbar ) 732 DO il=1,lml 733 DO ij=1,jml 734 DO ii=1,iml-1 735 champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij) 736 ENDDO 737 champ(iml,ij, il) = champ(1,ij, il) 738 ENDDO 739 ENDDO 740 CASE ('v') 741 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 742 CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2, 743 . lon_in2, lat_in2 , interbar ) 744 ENDIF 745 CALL start_inter_3d('V', iml, jml, lml, lon_in, 746 . lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar ) 747 DO il=1,lml 748 DO ij=1,jml 749 DO ii=1,iml-1 750 champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij) 751 ENDDO 752 champ(iml,ij, il) = champ(1,ij, il) 753 ENDDO 754 ENDDO 755 CASE ('t') 756 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 757 CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 , 758 . lon_in2, lat_in2 ,interbar ) 759 ENDIF 760 CALL start_inter_3d('TEMP', iml, jml, lml, lon_in, 761 . lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar ) 762 763 CASE ('tpot') 764 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 765 CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2 , 766 . lon_in2, lat_in2 , interbar ) 767 ENDIF 768 CALL start_inter_3d('TEMP', iml, jml, lml, lon_in, 769 . lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar ) 770 IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) ) 771 . THEN 772 DO il=1,lml 773 DO ij=1,jml 774 DO ii=1,iml-1 775 champ(ii,ij,il) = champ(ii,ij,il) * cpp 776 . / workvar(ii,ij,il) 777 ENDDO 778 champ(iml,ij,il) = champ(1,ij,il) 779 ENDDO 780 ENDDO 781 DO il=1,lml 782 xppn = SUM(aire(:,1)*champ(:,1,il))/apoln 783 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 784 champ(:,1,il) = xppn 785 champ(:,jml,il) = xpps 786 ENDDO 787 ELSE 788 WRITE(*,*)'Could not compute potential temperature as the' 789 WRITE(*,*)'Exner function is missing or constant.' 790 STOP 791 ENDIF 792 CASE ('q') 793 IF ( .NOT.ALLOCATED(psol_dyn)) THEN 794 CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 , 795 . lon_in2, lat_in2 , interbar ) 796 ENDIF 797 CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in, 798 . jml2, lon_in2, lat_in2, pls, champ, interbar ) 799 IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) ) 800 . THEN 801 DO il=1,lml 802 DO ij=1,jml 803 DO ii=1,iml-1 804 champ(ii,ij,il) = 0.01 * champ(ii,ij,il) * 805 . workvar(ii,ij,il) 806 ENDDO 807 champ(iml,ij,il) = champ(1,ij,il) 808 ENDDO 809 ENDDO 810 WHERE ( champ .LT. 0.) champ = 1.0E-10 811 DO il=1,lml 812 xppn = SUM(aire(:,1)*champ(:,1,il))/apoln 813 xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols 814 champ(:,1,il) = xppn 815 champ(:,jml,il) = xpps 816 ENDDO 817 ELSE 818 WRITE(*,*)'Could not compute specific humidity as the' 819 WRITE(*,*)'saturated humidity is missing or constant.' 820 STOP 821 ENDIF 822 CASE DEFAULT 823 WRITE(*,*) 'startget_dyn' 824 WRITE(*,*) 'No rule is present to extract variable ', 825 . varname(:LEN_TRIM(varname)),' from any data set' 826 STOP 827 END SELECT 828 ENDIF 829 END SUBROUTINE startget_dyn 830 ! 831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 832 ! 833 SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in,jml2,lon_in2 , 834 , lat_in2 , interbar ) 835 ! 836 INTEGER, INTENT(in) :: iml, jml, jml2 837 REAL, INTENT(in) :: lon_in(iml), lat_in(jml) 838 REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2) 839 LOGICAL interbar 840 ! 841 ! LOCAL 842 ! 843 REAL :: lev(1), date, dt 844 INTEGER :: itau(1) 845 INTEGER :: i, j 846 integer :: iret 847 ! 848 CHARACTER*120 :: physfname 849 LOGICAL :: check=.TRUE. 850 ! 851 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 852 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) 853 REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:) 854 REAL, ALLOCATABLE :: xppn(:), xpps(:) 855 LOGICAL :: allo 856 ! 857 ! 858 #include "dimensions.h" 859 #include "paramet.h" 860 #include "comgeom2.h" 861 862 CHARACTER*25 title 863 864 ! 865 physfname = 'ECDYN.nc' 866 ! 867 IF ( check ) WRITE(*,*) 'Opening the surface analysis' 868 ! 869 CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, 870 . ttm_dyn, fid_dyn) 871 IF ( check ) WRITE(*,*) 'Values read: ', iml_dyn, jml_dyn, 872 . llm_dyn, ttm_dyn 873 ! 874 ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret) 875 ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret) 876 ALLOCATE (levdyn_ini(llm_dyn), stat=iret) 877 ! 878 CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, 879 . lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, 880 . itau, date, dt, fid_dyn) 881 ! 882 883 allo = allocated (var_ana) 884 if (allo) then 885 DEALLOCATE(var_ana, stat=iret) 886 endif 887 ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret) 888 889 allo = allocated (lon_rad) 890 if (allo) then 891 DEALLOCATE(lon_rad, stat=iret) 892 endif 893 894 ALLOCATE(lon_rad(iml_dyn), stat=iret) 895 ALLOCATE(lon_ini(iml_dyn)) 896 897 IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 898 lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0 899 ELSE 900 lon_ini(:) = lon_dyn(:,1) 901 ENDIF 902 903 ALLOCATE(lat_rad(jml_dyn)) 904 ALLOCATE(lat_ini(jml_dyn)) 905 906 IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 907 lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0 908 ELSE 909 lat_ini(:) = lat_dyn(1,:) 910 ENDIF 911 ! 912 913 914 ALLOCATE(z(iml, jml)) 915 ALLOCATE(tmp_var(iml-1,jml)) 916 ! 917 CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn, 918 . 1, 1, var_ana) 919 c 920 title='Z' 921 CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini, 922 . lon_rad, lat_rad, var_ana, interbar ) 923 c 924 IF ( interbar ) THEN 925 WRITE(6,*) '-------------------------------------------------', 926 ,'--------------' 927 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 928 , ' pour Z $$$ ' 929 WRITE(6,*) '-------------------------------------------------', 930 ,'--------------' 931 CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad , 932 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var) 933 ELSE 934 CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana, 935 . iml-1, jml, lon_in, lat_in, tmp_var) 936 ENDIF 937 938 CALL gr_int_dyn(tmp_var, z, iml-1, jml) 939 ! 940 ALLOCATE(psol_dyn(iml, jml)) 941 ! 942 CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn, 943 . 1, 1, var_ana) 944 945 title='SP' 946 CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini, 947 . lon_rad, lat_rad, var_ana, interbar ) 948 949 IF ( interbar ) THEN 950 WRITE(6,*) '-------------------------------------------------', 951 ,'--------------' 952 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 953 , ' pour SP $$$ ' 954 WRITE(6,*) '-------------------------------------------------', 955 ,'--------------' 956 CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad , 957 , var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var) 958 ELSE 959 CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana, 960 . iml-1, jml, lon_in, lat_in, tmp_var ) 961 ENDIF 962 963 CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml) 964 ! 965 IF ( .NOT.ALLOCATED(tsol)) THEN 966 ! These variables may have been allocated by the need to 967 ! create a start field for them or by the varibale 968 ! coming out of the restart file. In case we dor have it we will initialize it. 969 ! 970 CALL start_init_phys( iml, jml, lon_in, lat_in,jml2,lon_in2, 971 . lat_in2 , interbar ) 972 ELSE 973 IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN 974 WRITE(*,*) 'start_init_dyn :' 975 WRITE(*,*) 'The temperature field we have does not ', 976 . 'have the right size' 977 STOP 978 ENDIF 979 ENDIF 980 IF ( .NOT.ALLOCATED(phis)) THEN 981 ! 982 ! These variables may have been allocated by the need to create a start field for them or by the varibale 983 ! coming out of the restart file. In case we dor have it we will initialize it. 984 ! 985 CALL start_init_orog( iml, jml, lon_in, lat_in, jml2, lon_in2 , 986 . lat_in2 , interbar ) 987 ! 988 ELSE 989 ! 990 IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN 991 ! 992 WRITE(*,*) 'start_init_dyn :' 993 WRITE(*,*) 'The orography field we have does not ', 994 . ' have the right size' 995 STOP 996 ENDIF 997 ! 998 ENDIF 999 ! 1000 ! PSOL is computed in Pascals 1001 ! 1002 ! 1003 DO j = 1, jml 1004 DO i = 1, iml-1 1005 psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j)) 1006 . /287.0/tsol(i,j)) 1007 ENDDO 1008 psol_dyn(iml,j) = psol_dyn(1,j) 1009 ENDDO 1010 ! 1011 ! 1012 ALLOCATE(xppn(iml-1)) 1013 ALLOCATE(xpps(iml-1)) 1014 ! 1015 DO i = 1, iml-1 1016 xppn(i) = aire( i,1) * psol_dyn( i,1) 1017 xpps(i) = aire( i,jml) * psol_dyn( i,jml) 1018 ENDDO 1019 ! 1020 DO i = 1, iml 1021 psol_dyn(i,1 ) = SUM(xppn)/apoln 1022 psol_dyn(i,jml) = SUM(xpps)/apols 1023 ENDDO 1024 ! 1025 RETURN 1026 ! 1027 END SUBROUTINE start_init_dyn 1028 ! 1029 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1030 ! 1031 SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, 1032 . lat_in, jml2, lon_in2, lat_in2, pls_in, var3d, interbar ) 1033 ! 1034 ! This subroutine gets a variables from a 3D file and does the interpolations needed 1035 ! 1036 ! 1037 ! ARGUMENTS 1038 ! 1039 CHARACTER*(*) :: varname 1040 INTEGER :: iml, jml, lml, jml2 1041 REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml) 1042 REAL :: lon_in2(iml) , lat_in2(jml2) 1043 REAL :: var3d(iml, jml, lml) 1044 LOGICAL interbar 1045 real chmin,chmax 1046 ! 1047 ! LOCAL 1048 ! 1049 CHARACTER*25 title 1050 INTEGER :: ii, ij, il, jsort,i,j,l 1051 REAL :: bx, by 1052 REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:) 1053 REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) , lev_dyn(:) 1054 REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:) 1055 REAL, ALLOCATABLE :: ax(:), ay(:), yder(:) 1056 ! REAL, ALLOCATABLE :: varrr(:,:,:) 1057 INTEGER, ALLOCATABLE :: lind(:) 1058 ! 1059 LOGICAL :: check = .TRUE. 1060 ! 1061 IF ( .NOT. ALLOCATED(var_ana3d)) THEN 1062 ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) 1063 ENDIF 1064 ! ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn)) 1065 ! 1066 ! 1067 IF ( check) WRITE(*,*) 'Going into flinget to extract the 3D ', 1068 . ' field.', fid_dyn 1069 IF ( check) WRITE(*,*) fid_dyn, varname, iml_dyn, jml_dyn, 1070 . llm_dyn,ttm_dyn 1071 ! 1072 CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, 1073 . ttm_dyn, 1, 1, var_ana3d) 1074 ! 1075 IF ( check) WRITE(*,*) 'Allocating space for the interpolation', 1076 . iml, jml, llm_dyn 1077 ! 1078 ALLOCATE(lon_rad(iml_dyn)) 1079 ALLOCATE(lon_ini(iml_dyn)) 1080 1081 IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 1082 lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0 1083 ELSE 1084 lon_ini(:) = lon_dyn(:,1) 1085 ENDIF 1086 1087 ALLOCATE(lat_rad(jml_dyn)) 1088 ALLOCATE(lat_ini(jml_dyn)) 1089 1090 ALLOCATE(lev_dyn(llm_dyn)) 1091 1092 IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN 1093 lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0 1094 ELSE 1095 lat_ini(:) = lat_dyn(1,:) 1096 ENDIF 1097 ! 1098 1099 CALL conf_dat3d ( varname,iml_dyn, jml_dyn, llm_dyn, lon_ini, 1100 . lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d , 1101 , interbar ) 1102 1103 ALLOCATE(var_tmp2d(iml-1, jml)) 1104 ALLOCATE(var_tmp3d(iml, jml, llm_dyn)) 1105 ALLOCATE(ax(llm_dyn)) 1106 ALLOCATE(ay(llm_dyn)) 1107 ALLOCATE(yder(llm_dyn)) 1108 ALLOCATE(lind(llm_dyn)) 1109 ! 1110 1111 DO il=1,llm_dyn 1112 ! 1113 IF( interbar ) THEN 1114 IF( il.EQ.1 ) THEN 1115 WRITE(6,*) '-------------------------------------------------', 1116 ,'--------------' 1117 WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ', 1118 , ' pour ', varname 1119 WRITE(6,*) '-------------------------------------------------', 1120 ,'--------------' 1121 ENDIF 1122 CALL inter_barxy ( iml_dyn, jml_dyn -1,lon_rad, lat_rad, 1123 , var_ana3d(:,:,il),iml-1, jml2, lon_in2, lat_in2,jml,var_tmp2d ) 1124 ELSE 1125 CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad, 1126 . var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d ) 1127 ENDIF 1128 ! 1129 CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml) 1130 ! 1131 ENDDO 1132 ! 1133 DO il=1,llm_dyn 1134 lind(il) = llm_dyn-il+1 1135 ENDDO 1136 ! 1137 c 1138 c ... Pour l'interpolation verticale ,on interpole du haut de l'atmosphere 1139 c vers le sol ... 1140 c 1141 DO ij=1,jml 1142 DO ii=1,iml-1 1143 ! 1144 ax(:) = lev_dyn(lind(:)) 1145 ay(:) = var_tmp3d(ii, ij, lind(:)) 1146 ! 1147 1148 CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder) 1149 ! 1150 DO il=1,lml 1151 bx = pls_in(ii, ij, il) 1152 CALL SPLINT(ax, ay, yder, llm_dyn, bx, by) 1153 var3d(ii, ij, il) = by 1154 ENDDO 1155 ! 1156 ENDDO 1157 var3d(iml, ij, :) = var3d(1, ij, :) 1158 ENDDO 1159 1160 do il=1,lml 1161 call minmax(iml*jml,var3d(1,1,il),chmin,chmax) 1162 SELECTCASE(varname) 1163 CASE('U') 1164 WRITE(*,*) ' U min max l ',il,chmin,chmax 1165 CASE('V') 1166 WRITE(*,*) ' V min max l ',il,chmin,chmax 1167 CASE('TEMP') 1168 WRITE(*,*) ' TEMP min max l ',il,chmin,chmax 1169 CASE('R') 1170 WRITE(*,*) ' R min max l ',il,chmin,chmax 1171 END SELECT 1172 enddo 1173 1174 DEALLOCATE(lon_rad) 1175 DEALLOCATE(lon_ini) 1176 DEALLOCATE(lat_rad) 1177 DEALLOCATE(lat_ini) 1178 DEALLOCATE(lev_dyn) 1179 DEALLOCATE(var_tmp2d) 1180 DEALLOCATE(var_tmp3d) 1181 DEALLOCATE(ax) 1182 DEALLOCATE(ay) 1183 DEALLOCATE(yder) 1184 DEALLOCATE(lind) 1185 1186 ! 1187 RETURN 1188 ! 1189 END SUBROUTINE start_inter_3d 1190 ! 530 CHARACTER(LEN=25) :: title 531 CHARACTER(LEN=120) :: physfname 532 LOGICAL :: check=.TRUE. 533 REAL :: date, dt 534 INTEGER, DIMENSION(1) :: itau 535 INTEGER :: i, j 536 REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana, z 537 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 538 REAL, DIMENSION(:), ALLOCATABLE :: xppn, xpps 539 !------------------------------------------------------------------------------- 540 541 !--- KINETIC ENERGY 542 physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0 543 IF(check) WRITE(lunout,*) 'Opening the surface analysis' 544 CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) 545 IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn 546 547 ALLOCATE(lat_dyn(iml_dyn,jml_dyn)) 548 ALLOCATE(lon_dyn(iml_dyn,jml_dyn)) 549 ALLOCATE(levdyn_ini(llm_dyn)) 550 CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,& 551 levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn) 552 553 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 554 ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) 555 lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad 556 lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad 557 558 ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) 559 560 !--- SURFACE GEOPOTENTIAL 561 title='Z' 562 ALLOCATE(z(iml,jml)) 563 CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana) 564 CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, & 565 var_ana, ibar) 566 CALL interp_startvar(title, ibar, .TRUE., & 567 iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 568 lon_in, lat_in, lon_in2, lat_in2, z) 569 570 !--- SURFACE PRESSURE 571 title='SP' 572 ALLOCATE(psol_dyn(iml,jml)) 573 CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana) 574 CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, & 575 var_ana, ibar) 576 CALL interp_startvar(title, ibar, .TRUE., & 577 iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & 578 lon_in, lat_in, lon_in2, lat_in2, psol_dyn) 579 580 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) 581 582 !--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE 583 IF(.NOT.ALLOCATED(tsol)) THEN 584 CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) 585 ELSE 586 IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN 587 WRITE(lunout,*)'start_init_dyn :' 588 WRITE(lunout,*)'The temperature field we have does not have the right size' 589 STOP 590 END IF 591 END IF 592 593 IF(.NOT.ALLOCATED(phis)) THEN 594 CALL start_init_orog(iml,jml,lon_in,lat_in) 595 ELSE 596 IF(SIZE(phis)/=SIZE(psol_dyn)) THEN 597 WRITE(lunout,*)'start_init_dyn :' 598 WRITE(lunout,*)'The orography field we have does not have the right size' 599 STOP 600 END IF 601 END IF 602 603 !--- PSOL IS COMPUTED IN PASCALS 604 DO j = 1, jml 605 DO i = 1, iml-1 606 psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j)) 607 END DO 608 psol_dyn(iml,j) = psol_dyn(1,j) 609 END DO 610 DEALLOCATE(z) 611 612 ALLOCATE(xppn(iml-1),xpps(iml-1)) 613 DO i = 1, iml-1 614 xppn(i) = aire( i,1) * psol_dyn( i,1) 615 xpps(i) = aire( i,jml) * psol_dyn( i,jml) 616 END DO 617 DO i = 1, iml 618 psol_dyn(i,1 ) = SUM(xppn)/apoln 619 psol_dyn(i,jml) = SUM(xpps)/apols 620 END DO 621 DEALLOCATE(xppn,xpps) 622 623 RETURN 624 625 END SUBROUTINE start_init_dyn 626 ! 627 !------------------------------------------------------------------------------- 628 629 630 !------------------------------------------------------------------------------- 631 ! 632 SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,& 633 lat_in2, pls_in, var3d, ibar) 634 ! 635 !------------------------------------------------------------------------------- 636 ! Arguments: 637 CHARACTER(LEN=*), INTENT(IN) :: varname 638 INTEGER, INTENT(IN) :: iml, jml, lml 639 REAL, DIMENSION(iml), INTENT(IN) :: lon_in 640 REAL, DIMENSION(jml), INTENT(IN) :: lat_in 641 INTEGER, INTENT(IN) :: jml2 642 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 643 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 644 REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls_in 645 REAL, DIMENSION(iml,jml,lml), INTENT(OUT) :: var3d 646 LOGICAL, INTENT(IN) :: ibar 647 !------------------------------------------------------------------------------- 648 ! Local variables: 649 #include "iniprint.h" 650 LOGICAL :: check=.TRUE. 651 REAL :: bx, by, chmin, chmax 652 INTEGER :: ii, ij, il 653 REAL, DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d 654 REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini 655 REAL, DIMENSION(:), ALLOCATABLE :: lev_dyn, ax, ay, yder 656 INTEGER, DIMENSION(:), ALLOCATABLE :: lind 657 !------------------------------------------------------------------------------- 658 IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D field.' 659 IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn 660 IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn 661 662 IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn)) 663 CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) 664 665 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 666 ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) 667 lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad 668 lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad 669 670 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 671 ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn)) 672 CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, & 673 levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar) 674 DEALLOCATE(lon_ini,lat_ini) 675 676 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro 677 ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) 678 DO il=1,llm_dyn 679 CALL interp_startvar(varname, ibar, il==1, & 680 iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2, & 681 lon_in, lat_in, lon_in2, lat_in2, var_tmp3d(:,:,il)) 682 END DO 683 DEALLOCATE(lon_rad,lat_rad) 684 685 ALLOCATE(lind(llm_dyn)) 686 DO il=1,llm_dyn 687 lind(il) = llm_dyn-il+1 688 END DO 689 690 !--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND 691 ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) 692 DO ij=1,jml 693 DO ii=1,iml-1 694 ax(:)=lev_dyn(lind(:)) 695 ay(:)=var_tmp3d(ii,ij,lind(:)) 696 CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder) 697 DO il=1,lml 698 bx=pls_in(ii,ij,il) 699 CALL SPLINT(ax, ay, yder, llm_dyn, bx, by) 700 var3d(ii,ij,il) = by 701 END DO 702 END DO 703 var3d(iml,ij,:) = var3d(1,ij,:) 704 END DO 705 DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind) 706 707 DO il=1,lml 708 CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax) 709 WRITE(lunout,*)' '//TRIM(varname)//' min max l ',il,chmin,chmax 710 END DO 711 712 RETURN 713 714 END SUBROUTINE start_inter_3d 715 ! 716 !------------------------------------------------------------------------------- 717 718 719 !------------------------------------------------------------------------------- 720 ! 721 SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj, lon, lat, vari, & 722 i1, j1, j2, lon1, lat1, lon2, lat2, varo) 723 ! 724 !------------------------------------------------------------------------------- 725 ! Arguments: 726 CHARACTER(LEN=*), INTENT(IN) :: vname 727 LOGICAL, INTENT(IN) :: ibar, ibeg 728 INTEGER, INTENT(IN) :: ii, jj 729 REAL, DIMENSION(ii), INTENT(IN) :: lon 730 REAL, DIMENSION(jj), INTENT(IN) :: lat 731 REAL, DIMENSION(ii,jj), INTENT(IN) :: vari 732 INTEGER, INTENT(IN) :: i1, j1, j2 733 REAL, DIMENSION(i1), INTENT(IN) :: lon1 734 REAL, DIMENSION(j1), INTENT(IN) :: lat1 735 REAL, DIMENSION(i1), INTENT(IN) :: lon2 736 REAL, DIMENSION(j2), INTENT(IN) :: lat2 737 REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo 738 !------------------------------------------------------------------------------- 739 ! Local variables: 740 #include "iniprint.h" 741 REAL, DIMENSION(i1-1,j1) :: vtmp 742 !------------------------------------------------------------------------------- 743 IF(ibar) THEN 744 IF(ibeg) THEN 745 WRITE(lunout,*) & 746 '---------------------------------------------------------------' 747 WRITE(lunout,*) & 748 '$$$ Utilisation de l interpolation barycentrique pour '//TRIM(vname)//' $$$' 749 WRITE(lunout,*) & 750 '---------------------------------------------------------------' 751 END IF 752 CALL inter_barxy(ii, jj-1, lon, lat, vari, i1-1, j2, lon2, lat2, j1, vtmp) 753 ELSE 754 CALL grille_m (ii, jj, lon, lat, vari, i1-1, j1, lon1, lat1, vtmp) 755 END IF 756 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 757 758 END SUBROUTINE interp_startvar 759 ! 760 !------------------------------------------------------------------------------- 761 1191 762 #endif 1192 763 ! of #ifdef CPP_EARTH 1193 END MODULE startvar 764 765 END MODULE startvar 766 ! 767 !*******************************************************************************
Note: See TracChangeset
for help on using the changeset viewer.
