source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/phyaqua_mod.F90 @ 3331

Last change on this file since 3331 was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 29.3 KB
Line 
1MODULE phyaqua_mod
2  ! Routines complementaires pour la physique planetaire.
3  IMPLICIT NONE
4
5CONTAINS
6
7  SUBROUTINE iniaqua(nlon, iflag_phys)
8
9    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10    ! Creation d'un etat initial et de conditions aux limites
11    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
12    ! du modele LMDZ dans sa version terrestre.
13    ! iflag_phys est un parametre qui controle
14    ! iflag_phys = N
15    ! de 100 a 199 : aqua planetes avec SST forcees
16    ! N-100 determine le type de SSTs
17    ! de 200 a 299 : terra planetes avec Ts calcule
18
19    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20
21    USE dimphy, ONLY: klon
22    USE geometry_mod, ONLY : latitude
23    USE surface_data, ONLY: type_ocean, ok_veget
24    USE pbl_surface_mod, ONLY: pbl_surface_init
25    USE fonte_neige_mod, ONLY: fonte_neige_init
26    USE phys_state_var_mod
27    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
28                                day_ini,day_end
29    USE indice_sol_mod
30    USE nrtype, ONLY: pi
31    USE ioipsl
32#ifdef ISO
33  USE infotrac_phy, ONLY: niso
34#endif
35    IMPLICIT NONE
36
37    include "YOMCST.h"
38    include "clesphys.h"
39    include "dimsoil.h"
40
41    INTEGER, INTENT (IN) :: nlon, iflag_phys
42    ! IM ajout latfi, lonfi
43!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
44
45    INTEGER type_profil, type_aqua
46
47    ! Ajouts initialisation des surfaces
48    REAL :: run_off_lic_0(nlon)
49    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
50    REAL :: tsoil(nlon, nsoilmx, nbsrf)
51    REAL :: tslab(nlon), seaice(nlon)
52    REAL fder(nlon)
53
54#ifdef ISO
55    REAL :: xtrun_off_lic_0(niso,nlon)
56    REAL :: xtsolsrf(niso,nlon, nbsrf), xtsnsrf(niso,nlon, nbsrf)
57    REAL :: Rland_ice(niso,nlon)
58#endif
59
60
61    ! Arguments :
62    ! -----------
63
64    ! integer radpas
65    INTEGER it, unit, i, k, itap
66
67    REAL airefi, zcufi, zcvfi
68
69    REAL rugos, albedo
70    REAL tsurf
71    REAL time, timestep, day, day0
72    REAL qsol_f
73    REAL rugsrel(nlon)
74    ! real zmea(nlon),zstd(nlon),zsig(nlon)
75    ! real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
76    ! real rlon(nlon),rlat(nlon)
77    LOGICAL alb_ocean
78    ! integer demih_pas
79
80    CHARACTER *80 ans, file_forctl, file_fordat, file_start
81    CHARACTER *100 file, var
82    CHARACTER *2 cnbl
83
84    REAL phy_nat(nlon, 360)
85    REAL phy_alb(nlon, 360)
86    REAL phy_sst(nlon, 360)
87    REAL phy_bil(nlon, 360)
88    REAL phy_rug(nlon, 360)
89    REAL phy_ice(nlon, 360)
90    REAL phy_fter(nlon, 360)
91    REAL phy_foce(nlon, 360)
92    REAL phy_fsic(nlon, 360)
93    REAL phy_flic(nlon, 360)
94
95    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
96
97    ! intermediate variables to use getin (need to be "save" to be shared by
98    ! all threads)
99    INTEGER, SAVE :: nbapp_rad_omp
100    REAL, SAVE :: co2_ppm_omp, solaire_omp
101    LOGICAL, SAVE :: alb_ocean_omp
102    REAL, SAVE :: rugos_omp
103    ! -------------------------------------------------------------------------
104    ! declaration pour l'appel a phyredem
105    ! -------------------------------------------------------------------------
106
107    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
108    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
109    ! real pbl_tke(nlon,llm,nbsrf)
110    ! real rain_fall(nlon),snow_fall(nlon)
111    ! real solsw(nlon), sollw(nlon),radsol(nlon)
112    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
113    ! real ratqs(nlon,llm)
114    ! real clwcon(nlon,llm)
115
116    INTEGER longcles
117    PARAMETER (longcles=20)
118    REAL clesphy0(longcles)
119
120
121    ! -----------------------------------------------------------------------
122    ! dynamial tendencies :
123    ! ---------------------
124
125    INTEGER l, ierr, aslun
126
127!    REAL longitude, latitude
128    REAL paire
129
130!    DATA latitude, longitude/48., 0./
131
132    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
133    ! INITIALISATIONS
134    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135
136    ! -----------------------------------------------------------------------
137    ! Initialisations  des constantes
138    ! -------------------------------
139
140
141    type_aqua = iflag_phys/100
142    type_profil = iflag_phys - type_aqua*100
143    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
144
145    IF (klon/=nlon) THEN
146      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
147      STOP 'probleme de dimensions dans iniaqua'
148    END IF
149    CALL phys_state_var_init(read_climoz)
150
151
152    read_climoz = 0
153    day0 = 217.
154    day = day0
155    it = 0
156    time = 0.
157
158    ! IM ajout latfi, lonfi
159!    rlatd = latfi
160!    rlond = lonfi
161!    rlat = rlatd*180./pi
162!    rlon = rlond*180./pi
163
164    ! -----------------------------------------------------------------------
165    ! initialisations de la physique
166    ! -----------------------------------------------------------------------
167
168    day_ini = day_ref
169    day_end = day_ini + ndays
170!    airefi = 1.
171!    zcufi = 1.
172!    zcvfi = 1.
173    !$OMP MASTER
174    nbapp_rad_omp = 24
175    CALL getin('nbapp_rad', nbapp_rad_omp)
176    !$OMP END MASTER
177    !$OMP BARRIER
178    nbapp_rad = nbapp_rad_omp
179
180    ! ---------------------------------------------------------------------
181    ! Creation des conditions aux limites:
182    ! ------------------------------------
183    ! Initialisations des constantes
184    ! Ajouter les manquants dans planete.def... (albedo etc)
185    !$OMP MASTER
186    co2_ppm_omp = 348.
187    CALL getin('co2_ppm', co2_ppm_omp)
188    solaire_omp = 1365.
189    CALL getin('solaire', solaire_omp)
190    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
191    ! type_aqua
192    alb_ocean_omp = .TRUE.
193    CALL getin('alb_ocean', alb_ocean_omp)
194    !$OMP END MASTER
195    !$OMP BARRIER
196    co2_ppm = co2_ppm_omp
197    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
198    solaire = solaire_omp
199    WRITE (*, *) 'iniaqua: solaire=', solaire
200    alb_ocean = alb_ocean_omp
201    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
202
203    radsol = 0.
204    qsol_f = 10.
205
206    ! Conditions aux limites:
207    ! -----------------------
208
209    qsol(:) = qsol_f
210    rugsrel = 0.0 ! (rugsrel = rugoro)
211    rugoro = 0.0
212    u_ancien = 0.0
213    v_ancien = 0.0
214    agesno = 50.0
215    ! Relief plat
216    zmea = 0.
217    zstd = 0.
218    zsig = 0.
219    zgam = 0.
220    zthe = 0.
221    zpic = 0.
222    zval = 0.
223
224    ! Une seule surface
225    pctsrf = 0.
226    IF (type_aqua==1) THEN
227      rugos = 1.E-4
228      albedo = 0.19
229      pctsrf(:, is_oce) = 1.
230    ELSE IF (type_aqua==2) THEN
231      rugos = 0.03
232      albedo = 0.1
233      pctsrf(:, is_ter) = 1.
234    END IF
235
236    !$OMP MASTER
237    rugos_omp = rugos
238    CALL getin('rugos', rugos_omp)
239    !$OMP END MASTER
240    !$OMP BARRIER
241    rugos = rugos_omp
242    WRITE (*, *) 'iniaqua: rugos=', rugos
243    zmasq(:) = pctsrf(:, is_ter)
244
245    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
246    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
247
248    ! Si alb_ocean on calcule un albedo oceanique moyen
249    ! if (alb_ocean) then
250    ! Voir pourquoi on avait ca.
251    ! CALL ini_alb_oce(phy_alb)
252    ! else
253    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
254    ! endif !alb_ocean
255
256    DO i = 1, 360
257      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
258      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
259
260      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
261      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
262      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
263      phy_ice(:, i) = 0.0 ! fraction de glace (?)
264      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
265      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
266      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
267      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
268    END DO
269    ! IM calcul profil sst
270    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
271
272    CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
273      phy_fter, phy_foce, phy_flic, phy_fsic)
274
275
276    ! ---------------------------------------------------------------------
277    ! Ecriture de l'etat initial:
278    ! ---------------------------
279
280
281    ! Ecriture etat initial physique
282
283    timestep = pdtphys
284    radpas = nint(rday/timestep/float(nbapp_rad))
285
286    DO i = 1, longcles
287      clesphy0(i) = 0.
288    END DO
289    clesphy0(1) = float(iflag_con)
290    clesphy0(2) = float(nbapp_rad)
291    ! IF( cycle_diurne  ) clesphy0(3) =  1.
292    clesphy0(3) = 1. ! cycle_diurne
293    clesphy0(4) = 1. ! soil_model
294    clesphy0(5) = 1. ! new_oliq
295    clesphy0(6) = 0. ! ok_orodr
296    clesphy0(7) = 0. ! ok_orolf
297    clesphy0(8) = 0. ! ok_limitvrai
298
299
300    ! =======================================================================
301    ! Profils initiaux
302    ! =======================================================================
303
304    ! On initialise les temperatures de surfaces comme les sst
305    DO i = 1, nlon
306      ftsol(i, :) = phy_sst(i, 1)
307      tsoil(i, :, :) = phy_sst(i, 1)
308      tslab(i) = phy_sst(i, 1)
309    END DO
310
311    falbe(:, :) = albedo
312    falblw(:, :) = albedo
313    rain_fall(:) = 0.
314    snow_fall(:) = 0.
315    solsw(:) = 0.
316    sollw(:) = 0.
317    radsol(:) = 0.
318
319    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
320    ! intialisation bidon mais pas grave
321    t_ancien(:, :) = 0.
322    q_ancien(:, :) = 0.
323    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324    rnebcon = 0.
325    ratqs = 0.
326    clwcon = 0.
327    pbl_tke = 1.E-8
328
329    ! variables supplementaires pour appel a plb_surface_init
330    fder(:) = 0.
331    seaice(:) = 0.
332    run_off_lic_0 = 0.
333    fevap = 0.
334#ifdef ISO
335    xtrun_off_lic_0 = 0.
336#endif
337
338    ! Initialisations necessaires avant phyredem
339    type_ocean = 'force'
340    CALL fonte_neige_init(run_off_lic_0 &
341#ifdef ISO
342    &  ,xtrun_off_lic_0 &
343#endif     
344    &   )
345    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
346    snsrf(:, :) = 0. ! couverture de neige des sous surface
347    z0m(:, :) = rugos ! couverture de neige des sous surface
348    z0h=z0m
349#ifdef ISO
350    xtsolsrf(:, :,:) = qsol(1)
351    xtsnsrf(:, :,:) = 0.   
352#endif
353
354    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil &
355#ifdef ISO       
356     &   ,xtsnsrf,Rland_ice &
357#endif
358     &  )
359
360    PRINT *, 'iniaqua: before phyredem'
361
362    falb1 = albedo
363    falb2 = albedo
364    zmax0 = 0.
365    f0 = 0.
366    sig1 = 0.
367    w01 = 0.
368    wake_deltat = 0.
369    wake_deltaq = 0.
370    wake_s = 0.
371    wake_cstar = 0.
372    wake_pe = 0.
373    wake_fip = 0.
374    fm_therm = 0.
375    entr_therm = 0.
376    detr_therm = 0.
377
378
379    CALL phyredem('startphy.nc')
380
381    PRINT *, 'iniaqua: after phyredem'
382    CALL phys_state_var_end
383
384    RETURN
385  END SUBROUTINE iniaqua
386
387
388  ! ====================================================================
389  ! ====================================================================
390  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
391    USE dimphy
392    IMPLICIT NONE
393    ! ====================================================================
394    ! =============================================================
395    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
396    ! Auteur : A. Campoy et F. Hourdin
397    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
398    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
399    ! connaissant la declinaison, la latitude et la longitude.
400
401    ! Dans cette version particuliere, on calcule le rayonnement
402    ! moyen sur l'année à chaque latitude.
403    ! angle zenithal calculé pour obtenir un
404    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
405    ! en moyenne annuelle.
406    ! Spécifique de la terre. Utilisé pour les aqua planetes.
407
408    ! Rque   : Different de la routine angle en ce sens que zenang
409    ! fournit des moyennes de pmu0 et non des valeurs
410    ! instantanees, du coup frac prend toutes les valeurs
411    ! entre 0 et 1.
412    ! Date   : premiere version le 13 decembre 1994
413    ! revu pour  GCM  le 30 septembre 1996
414    ! ===============================================================
415    ! longi----INPUT : la longitude vraie de la terre dans son plan
416    ! solaire a partir de l'equinoxe de printemps (degre)
417    ! gmtime---INPUT : temps universel en fraction de jour
418    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
419    ! lat------INPUT : latitude en degres
420    ! long-----INPUT : longitude en degres
421    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
422    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
423    ! ================================================================
424    include "YOMCST.h"
425    ! ================================================================
426    LOGICAL cycle_diurne
427    REAL gmtime
428    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
429    ! ================================================================
430    INTEGER i
431    REAL gmtime1, gmtime2
432    REAL pi_local
433
434
435    REAL rmu0m(klon), rmu0a(klon)
436
437
438    pi_local = 4.0*atan(1.0)
439
440    ! ================================================================
441    ! Calcul de l'angle zenithal moyen sur la journee
442    ! ================================================================
443
444    DO i = 1, klon
445      fract(i) = 1.
446      ! Calcule du flux moyen
447      IF (abs(rlat(i))<=28.75) THEN
448        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
449      ELSE IF (abs(rlat(i))<=43.75) THEN
450        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
451      ELSE IF (abs(rlat(i))<=71.25) THEN
452        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
453      ELSE
454        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
455      END IF
456    END DO
457
458    ! ================================================================
459    ! Avec ou sans cycle diurne
460    ! ================================================================
461
462    IF (cycle_diurne) THEN
463
464      ! On redecompose flux  au sommet suivant un cycle diurne idealise
465      ! identique a toutes les latitudes.
466
467      DO i = 1, klon
468        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
469        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
470          rmu0a(i)/sqrt(2.)
471      END DO
472
473      DO i = 1, klon
474        IF (rmu0(i)<=0.) THEN
475          rmu0(i) = 0.
476          fract(i) = 0.
477        ELSE
478          fract(i) = 1.
479        END IF
480      END DO
481
482      ! Affichage de l'angel zenitale
483      ! print*,'************************************'
484      ! print*,'************************************'
485      ! print*,'************************************'
486      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
487      ! print*,'rmu0m=',rmu0m(i)
488      ! print*,'rmu0a=',rmu0a(i)
489      ! print*,'rmu0=',rmu0(i)
490
491    ELSE
492
493      DO i = 1, klon
494        fract(i) = 0.5
495        rmu0(i) = rmu0m(i)*2.
496      END DO
497
498    END IF
499
500    RETURN
501  END SUBROUTINE zenang_an
502
503  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
504
505  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
506      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
507
508    USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_omp_root
509    USE mod_grid_phy_lmdz, ONLY: klon_glo
510    USE mod_phys_lmdz_transfert_para, ONLY: gather
511    IMPLICIT NONE
512    include "netcdf.inc"
513
514    INTEGER, INTENT (IN) :: klon
515    REAL, INTENT (IN) :: phy_nat(klon, 360)
516    REAL, INTENT (IN) :: phy_alb(klon, 360)
517    REAL, INTENT (IN) :: phy_sst(klon, 360)
518    REAL, INTENT (IN) :: phy_bil(klon, 360)
519    REAL, INTENT (IN) :: phy_rug(klon, 360)
520    REAL, INTENT (IN) :: phy_ice(klon, 360)
521    REAL, INTENT (IN) :: phy_fter(klon, 360)
522    REAL, INTENT (IN) :: phy_foce(klon, 360)
523    REAL, INTENT (IN) :: phy_flic(klon, 360)
524    REAL, INTENT (IN) :: phy_fsic(klon, 360)
525
526    REAL :: phy_glo(klon_glo, 360) ! temporary variable, to store phy_***(:)
527      ! on the whole physics grid
528    INTEGER :: k
529    INTEGER ierr
530    INTEGER dimfirst(3)
531    INTEGER dimlast(3)
532
533    INTEGER nid, ndim, ntim
534    INTEGER dims(2), debut(2), epais(2)
535    INTEGER id_tim
536    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
537    INTEGER id_fter, id_foce, id_fsic, id_flic
538
539    IF (is_mpi_root .AND. is_omp_root) THEN
540
541      PRINT *, 'writelim: Ecriture du fichier limit'
542
543      ierr = nf_create('limit.nc', nf_clobber, nid)
544
545      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
546        'Fichier conditions aux limites')
547      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
548      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
549      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
550
551      dims(1) = ndim
552      dims(2) = ntim
553
554#ifdef NC_DOUBLE
555      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
556#else
557      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
558#endif
559      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
560
561#ifdef NC_DOUBLE
562      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
563#else
564      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
565#endif
566      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
567        'Nature du sol (0,1,2,3)')
568
569#ifdef NC_DOUBLE
570      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
571#else
572      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
573#endif
574      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
575        'Temperature superficielle de la mer')
576
577#ifdef NC_DOUBLE
578      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
579#else
580      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
581#endif
582      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
583        'Reference flux de chaleur au sol')
584
585#ifdef NC_DOUBLE
586      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
587#else
588      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
589#endif
590      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
591
592#ifdef NC_DOUBLE
593      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
594#else
595      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
596#endif
597      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
598
599#ifdef NC_DOUBLE
600      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
601#else
602      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
603#endif
604      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
605#ifdef NC_DOUBLE
606      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
607#else
608      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
609#endif
610      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
611#ifdef NC_DOUBLE
612      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
613#else
614      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
615#endif
616      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
617#ifdef NC_DOUBLE
618      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
619#else
620      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
621#endif
622      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
623
624      ierr = nf_enddef(nid)
625      IF (ierr/=nf_noerr) THEN
626        WRITE (*, *) 'writelim error: failed to end define mode'
627        WRITE (*, *) nf_strerror(ierr)
628      END IF
629
630
631      ! write the 'times'
632      DO k = 1, 360
633#ifdef NC_DOUBLE
634        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
635#else
636        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
637#endif
638        IF (ierr/=nf_noerr) THEN
639          WRITE (*, *) 'writelim error with temps(k),k=', k
640          WRITE (*, *) nf_strerror(ierr)
641        END IF
642      END DO
643
644    END IF ! of if (is_mpi_root.and.is_omp_root)
645
646    ! write the fields, after having collected them on master
647
648    CALL gather(phy_nat, phy_glo)
649    IF (is_mpi_root .AND. is_omp_root) THEN
650#ifdef NC_DOUBLE
651      ierr = nf_put_var_double(nid, id_nat, phy_glo)
652#else
653      ierr = nf_put_var_real(nid, id_nat, phy_glo)
654#endif
655      IF (ierr/=nf_noerr) THEN
656        WRITE (*, *) 'writelim error with phy_nat'
657        WRITE (*, *) nf_strerror(ierr)
658      END IF
659    END IF
660
661    CALL gather(phy_sst, phy_glo)
662    IF (is_mpi_root .AND. is_omp_root) THEN
663#ifdef NC_DOUBLE
664      ierr = nf_put_var_double(nid, id_sst, phy_glo)
665#else
666      ierr = nf_put_var_real(nid, id_sst, phy_glo)
667#endif
668      IF (ierr/=nf_noerr) THEN
669        WRITE (*, *) 'writelim error with phy_sst'
670        WRITE (*, *) nf_strerror(ierr)
671      END IF
672    END IF
673
674    CALL gather(phy_bil, phy_glo)
675    IF (is_mpi_root .AND. is_omp_root) THEN
676#ifdef NC_DOUBLE
677      ierr = nf_put_var_double(nid, id_bils, phy_glo)
678#else
679      ierr = nf_put_var_real(nid, id_bils, phy_glo)
680#endif
681      IF (ierr/=nf_noerr) THEN
682        WRITE (*, *) 'writelim error with phy_bil'
683        WRITE (*, *) nf_strerror(ierr)
684      END IF
685    END IF
686
687    CALL gather(phy_alb, phy_glo)
688    IF (is_mpi_root .AND. is_omp_root) THEN
689#ifdef NC_DOUBLE
690      ierr = nf_put_var_double(nid, id_alb, phy_glo)
691#else
692      ierr = nf_put_var_real(nid, id_alb, phy_glo)
693#endif
694      IF (ierr/=nf_noerr) THEN
695        WRITE (*, *) 'writelim error with phy_alb'
696        WRITE (*, *) nf_strerror(ierr)
697      END IF
698    END IF
699
700    CALL gather(phy_rug, phy_glo)
701    IF (is_mpi_root .AND. is_omp_root) THEN
702#ifdef NC_DOUBLE
703      ierr = nf_put_var_double(nid, id_rug, phy_glo)
704#else
705      ierr = nf_put_var_real(nid, id_rug, phy_glo)
706#endif
707      IF (ierr/=nf_noerr) THEN
708        WRITE (*, *) 'writelim error with phy_rug'
709        WRITE (*, *) nf_strerror(ierr)
710      END IF
711    END IF
712
713    CALL gather(phy_fter, phy_glo)
714    IF (is_mpi_root .AND. is_omp_root) THEN
715#ifdef NC_DOUBLE
716      ierr = nf_put_var_double(nid, id_fter, phy_glo)
717#else
718      ierr = nf_put_var_real(nid, id_fter, phy_glo)
719#endif
720      IF (ierr/=nf_noerr) THEN
721        WRITE (*, *) 'writelim error with phy_fter'
722        WRITE (*, *) nf_strerror(ierr)
723      END IF
724    END IF
725
726    CALL gather(phy_foce, phy_glo)
727    IF (is_mpi_root .AND. is_omp_root) THEN
728#ifdef NC_DOUBLE
729      ierr = nf_put_var_double(nid, id_foce, phy_glo)
730#else
731      ierr = nf_put_var_real(nid, id_foce, phy_glo)
732#endif
733      IF (ierr/=nf_noerr) THEN
734        WRITE (*, *) 'writelim error with phy_foce'
735        WRITE (*, *) nf_strerror(ierr)
736      END IF
737    END IF
738
739    CALL gather(phy_fsic, phy_glo)
740    IF (is_mpi_root .AND. is_omp_root) THEN
741#ifdef NC_DOUBLE
742      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
743#else
744      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
745#endif
746      IF (ierr/=nf_noerr) THEN
747        WRITE (*, *) 'writelim error with phy_fsic'
748        WRITE (*, *) nf_strerror(ierr)
749      END IF
750    END IF
751
752    CALL gather(phy_flic, phy_glo)
753    IF (is_mpi_root .AND. is_omp_root) THEN
754#ifdef NC_DOUBLE
755      ierr = nf_put_var_double(nid, id_flic, phy_glo)
756#else
757      ierr = nf_put_var_real(nid, id_flic, phy_glo)
758#endif
759      IF (ierr/=nf_noerr) THEN
760        WRITE (*, *) 'writelim error with phy_flic'
761        WRITE (*, *) nf_strerror(ierr)
762      END IF
763    END IF
764
765    ! close file:
766    IF (is_mpi_root .AND. is_omp_root) THEN
767      ierr = nf_close(nid)
768    END IF
769
770  END SUBROUTINE writelim
771
772  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
773
774  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
775    USE dimphy
776    IMPLICIT NONE
777
778    INTEGER nlon, type_profil, i, k, j
779    REAL :: rlatd(nlon), phy_sst(nlon, 360)
780    INTEGER imn, imx, amn, amx, kmn, kmx
781    INTEGER p, pplus, nlat_max
782    PARAMETER (nlat_max=72)
783    REAL x_anom_sst(nlat_max)
784
785    IF (klon/=nlon) STOP 'probleme de dimensions dans iniaqua'
786    WRITE (*, *) ' profil_sst: type_profil=', type_profil
787    DO i = 1, 360
788      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
789
790      ! Rajout fbrlmd
791
792      IF (type_profil==1) THEN
793        ! Méthode 1 "Control" faible plateau à l'Equateur
794        DO j = 1, klon
795          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
796          ! PI/3=1.047197551
797          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
798            phy_sst(j, i) = 273.
799          END IF
800        END DO
801      END IF
802      IF (type_profil==2) THEN
803        ! Méthode 2 "Flat" fort plateau à l'Equateur
804        DO j = 1, klon
805          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
806          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
807            phy_sst(j, i) = 273.
808          END IF
809        END DO
810      END IF
811
812
813      IF (type_profil==3) THEN
814        ! Méthode 3 "Qobs" plateau réel à l'Equateur
815        DO j = 1, klon
816          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
817            rlatd(j))**4)
818          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
819            phy_sst(j, i) = 273.
820          END IF
821        END DO
822      END IF
823
824      IF (type_profil==4) THEN
825        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
826        DO j = 1, klon
827          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
828            rlatd(j))**4)
829          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
830            phy_sst(j, i) = 273.
831          END IF
832        END DO
833      END IF
834
835      IF (type_profil==5) THEN
836        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
837        DO j = 1, klon
838          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
839            *rlatd(j))**4)
840          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
841            phy_sst(j, i) = 273. + 2.
842          END IF
843
844        END DO
845      END IF
846
847      IF (type_profil==6) THEN
848        ! Méthode 6 "cst" valeur constante de SST
849        DO j = 1, klon
850          phy_sst(j, i) = 288.
851        END DO
852      END IF
853
854
855      IF (type_profil==7) THEN
856        ! Méthode 7 "cst" valeur constante de SST +2
857        DO j = 1, klon
858          phy_sst(j, i) = 288. + 2.
859        END DO
860      END IF
861
862      p = 0
863      IF (type_profil==8) THEN
864        ! Méthode 8 profil anomalies SST du modèle couplé AR4
865        DO j = 1, klon
866          IF (rlatd(j)==rlatd(j-1)) THEN
867            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
868              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
869            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
870              phy_sst(j, i) = 273. + x_anom_sst(pplus)
871            END IF
872          ELSE
873            p = p + 1
874            pplus = 73 - p
875            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
876              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
877            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
878              phy_sst(j, i) = 273. + x_anom_sst(pplus)
879            END IF
880            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
881          END IF
882        END DO
883      END IF
884
885      IF (type_profil==9) THEN
886        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
887        DO j = 1, klon
888          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
889            *rlatd(j))**4)
890          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
891            phy_sst(j, i) = 273. - 2.
892          END IF
893        END DO
894      END IF
895
896
897      IF (type_profil==10) THEN
898        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
899        DO j = 1, klon
900          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
901            *rlatd(j))**4)
902          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
903            phy_sst(j, i) = 273. + 4.
904          END IF
905        END DO
906      END IF
907
908      IF (type_profil==11) THEN
909        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
910        DO j = 1, klon
911          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
912            rlatd(j))**4)
913          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
914            phy_sst(j, i) = 273.
915          END IF
916        END DO
917      END IF
918
919      IF (type_profil==12) THEN
920        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
921        DO j = 1, klon
922          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
923            *rlatd(j))**4)
924          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
925            phy_sst(j, i) = 273. + 4.
926          END IF
927        END DO
928      END IF
929
930      IF (type_profil==13) THEN
931        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
932        DO j = 1, klon
933          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
934            rlatd(j))**4)
935          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
936            phy_sst(j, i) = 273.
937          END IF
938        END DO
939      END IF
940
941      IF (type_profil==14) THEN
942        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
943        DO j = 1, klon
944          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
945            *rlatd(j))**4)
946          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
947            phy_sst(j, i) = 273.
948          END IF
949        END DO
950      END IF
951
952      if (type_profil.EQ.20) then
953      print*,'Profile SST 20'
954!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
955
956      do j=1,klon
957        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
958      enddo
959      endif
960
961      if (type_profil.EQ.21) then
962      print*,'Profile SST 21'
963!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
964      do j=1,klon
965        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
966      enddo
967      endif
968
969
970
971    END DO
972
973    ! IM beg : verif profil SST: phy_sst
974    amn = min(phy_sst(1,1), 1000.)
975    amx = max(phy_sst(1,1), -1000.)
976    imn = 1
977    kmn = 1
978    imx = 1
979    kmx = 1
980    DO k = 1, 360
981      DO i = 2, nlon
982        IF (phy_sst(i,k)<amn) THEN
983          amn = phy_sst(i, k)
984          imn = i
985          kmn = k
986        END IF
987        IF (phy_sst(i,k)>amx) THEN
988          amx = phy_sst(i, k)
989          imx = i
990          kmx = k
991        END IF
992      END DO
993    END DO
994
995    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
996    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
997    ! IM end : verif profil SST: phy_sst
998
999    RETURN
1000  END SUBROUTINE profil_sst
1001
1002END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.